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CHAPTER  I 


INTRODUCTION 


1.1  Fourth  State  of  Matter. 

Solids,  liquids,  and  gases  are  the  well  known  three 
states  of  matter.  "Solids"  represent  the  state  in  which  the 
relative  motion  of  molecules  is  restricted  and  tends  to 
retain  a  definite  fixed  position  relative  to  each  other, 
which  results  in  a  crystalline  structure1.  'Liquids'' 
represent  the  state  in  which  molecules  are  relatively  free 
to  change  their  positions  about  one  another,  while 
restricted  by  cohesive  forces  such  that  they  maintain  a 
relatively  fixed  volume1.  And,  "gases"  represent  the  state 
in  whicn  molecules  are  virtually  unrestricted  by  cohesive 
forces,  resuiting  in  neither  definite  shape  nor  volume1. 

A  little  over  a  hundred  years  ago,  it  was  believed  that 
another  state  of  matter  was  observed.  This  state  was 
characterized  by  an  enclosed  electrically  neutral  collection 
of  ions,  electrons,  neutral  atoms  and  molecules.  It  was 
further  characterized  by  relatively  large  intermolecular 
distances  and  large  internal  energy  in  the  particles1.  The 
cnaracteristics  of  this  newly  discovered  matter  do  not  fit 
any  of  the  known  three  states  of  matter  provided  that  gas 
does  not  represent  a  state  of  a  collection  of  neutrals, 
ions,  and  electrons.  Thus,  it  is  clear  that  matter  exists 
in  more  than  three  states.  Commonly  referred  to  as  the 
"fourth  state  of  matter,"  this  new  state  is  called  the 
plasma  state. 

Plasmas  exist  in  many  forms.  These  forms  range  from 
the  hot  classical  plasma  found  in  the  magnetospheres  of 


1 


2 


pulsars  to  the  cold,  dense  degenerate  quantum  electron 

2 

plasma  of  a  white  dwarf  .  Usually,  plasmas  are 
characterized  by  a  high  degree  of  electrical  conductivity3. 

A  typical  gas  acts  as  a  good  electrical  insulator.  If 
one  were  to  produce  a  plasma  out  of  a  normal  gas,  one  would 
nave  a  useful  electrical  component  -  a  good  conductor. 

With  this  potential  usefulness,  plasmas  have  been 
artificially  produced  in  the  laboratory.  Some  of  the  ways 
in  which  a  plasma  may  be  produced  are  shock,  spark 
discharge,  arc  discharge,  nuclear  reaction,  chemical 
reaction  (of  high  specific  energy),  and  bombardment  by 

3 

electromagnetic  fields  . 

1.2  Plasma  Applications. 

Plasmas  can  be  used  for  many  applications.  Some  of 
these  include  production  of  nuclear  fuels,  research  and 

diagnostics  of  medicine,  research  in  agriculture,  and 

2 

environmental  tracking  of  pollutants  .  Compared  to 
conventional  metal  combinations,  the  use  of  plasma 
thermocouples  allows  us  to  extract  more  thermoelectric 

3 

power  from  nuclear  reactors  . 

For  military  purposes,  plasmas  can  be  used  for 

...  .  ,  .  2 3 

nitration  systems  in  a  toxic  chemical  environment  .  Also, 

plasmas  provide  useful  emission  sources  for  producing 

emission  spectra  for  chemical  analysis.  Additionally, 

plasmas  provide  a  novel  and  useful  role  in  jet  propulsion 

2 

for  space  flight  . 

1.3  Rocket  Propulsion. 

There  are  three  major  types  of  rocket  thrusters: 
chemical,  nuclear,  and  electrical. 
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chemical  rocket  thrusters  are  the  most  commonly  used 
type  of  tnruster.  Through  chemical  reactions  (combustion), 
energy  of  a  liquid  propellant  is  transformed  into  internal 
energy  in  the  form  of  hot,  high  pressure  gas.  These  gases 
expand  through  a  nozzle  and  form  a  propulsive  jet* . 

Nuclear  rocket  thrusters  require  that  a  working  fluid 
be  present.  Energy  is  transferred  to  the  working  fluid 
through  temperature  increases  caused  by  the  nuclear  energy. 
This  worKing  fluid  can  perform  one  of  two  different 
functions.  First,  it  can  expand  through  a  nozzle  (similar 
to  chemical  propulsion)  and  act  as  a  propulsive  jet. 

Second,  it  can  be  converted  to  electrical  energy  lor  power 
in  an  electrical  propulsion  process.  Unfortunately,  the 
practical  and  political  difficulties  involved  in  using 
nuciear  rocket  propulsion  are  likely  to  hinder  its 
development  to  widespread  use*. 

Electrical  rocket  thrusters  differ  substantially  from 

the  other  two  in  that  the  thruster  must  include  some  type 

of  power  plant  to  produce  the  electrical  power  for  the 

thrust  unit.  This  type  of  propulsion  is  useful  for 

steering  spacecraft  already  in  space.  However,  electrical 

thrusters  are  not  useful  where  gravity  has  a  significant 
5 

counteraction  . 

1.4  Microwave  Induction. 

Microwave  induced  plasmas  are  very  efficient  for  uses 
in  jet  propulsion.  Production  of  these  piasraas  involves 
using  plasma  columns  the  size  of  conventional  resonance 
cavities.  These  cavities  are  stable,  reproducible,  and 
quiescent.  These  plasmas  develop  as  a  result  of  surface 
wave  propagation  and  are  characterized  by  ion  immobility. 
The  major  physical  processes  governing  the  discharge  are: 
(a)  discharge  conditions  (such  as  the  nature  of  the  gas  i  , 


4 

(b)  gas  pressure,  (c)  dimensions  and  material  of  the  vessel, 
(d)  frequency  of  the  electromagnetic  field,  and  (e)  the 
power  transferred  to  the  plasma7. 

For  microwave  plasma  electrothermal  rocket  thrusters, 
pressures  of  about  one  atmosphere  and  temperatures  of  about 
1500  degrees  Kelvin  are  being  investigated.  In  the 
electromagnetic  environment,  free  electrons  are  accelerated 
about  the  heavier  and  neutral  molecules.  These  electrons 
coll ide  with  other  elements  or  molecules  of  the  gas  to 
cause  them  to  ionize  as  previously  bound  electrons  are 
stripped  off.  In  essence,  kinetic  energy  is  transferred 
from  trie  accelerated  electrons  to  the  gas. 

Figure  1.1  illustrates  the  various  discharge 
properties  within  a  microwave  system.  The  cold  propellant 
receives  microwave  energy  resulting  in  the  production  of  a 
plasma.  This  plasma  gives  off  radiation  and  heat  energy. 

The  excited  species  flows  away  from  the  plasma  while  the 
cold  species  flows  towards  it.  Finally,  the  plasma  excited 
propellant  is  recombined  downstream  of  the  plasma  with 
increased  kinetic  energy.  This  thermalized  propellant 
exits  through  a  nozzle  as  propulsion  thrust. 

1.5  Research  Direction. 

The  research  in  this  thesis  is  divided  into  three 
sections.  The  first  section  covers  a  theoretical 
background  in  the  microscopic  analysis  of  plasma  physics. 
Additionally,  a  brief  discussion  is  provided  for  potential 
pxasma  electric  propulsion. 

The  next  section  covers  an  experimental  study  of  a 
microwave  induced  plasma  in  a  cylindrical  cavity.  The 
macroscopic  power  distribution,  plasma  dimensions,  and 
electron  temperatures  are  determined  from  measurements  taken 


Figure  1.1  Microwave  Plasma  Discharge  Properties 
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of  helium  and  nitrogen  plasmas.  The  pressure  range  is 
between  200-1000  torr  with  gas  flow  rates  from  0-2000  SCCM. 
The  input  power  varies  between  200-275  watts. 

The  last  section  contains  a  computer  analysis  of  the 
electron  diffusion  within  the  plasma  area.  A  computer 
program  was  written  to  calculate  the  normalized  electron 
density  using  the  equations  derived  in  the  theoretical 
section.  The  computer  simulation  was  used  to  calculate  the 
electron  gradient  within  the  plasma  discharge  region.  The 
major  assumptions  for  this  model  are  constant  temperature 
and  equal  ion  4  electron  densities. 


CHAPTER  II 


THEORY 


2.1  Introduction. 

Modeling  transport  properties  of  the  plasma  region  in 
the  microwave  electrothermal  thruster  is  a  dynamic  process 
The  first  section  of  this  chapter  covers  rocket  propulsion 
physics  and  identifies  important  concepts.  The  next 
section  deals  with  the  complexities  involved  in  attempting 
to  characterize  the  plasma  region  of  the  proposed  thruster 
system.  The  latter  section  covers  the  theory  involved  in 
the  experimental  system. 

2.2  Rocket  Propulsion  Physics. 

Two  different  aspects,  thrust  and  geo-centric 
trajectory  (orbit),  will  be  discussed.  Not  much  detail 
will  be  given;  but,  several  terms  will  be  defined. 

2.2.1  Thrust . 


Thrust  is  defined  as  the  force  required  to  maintain 

velocity  of  expelling  propellents  as  mass  changes  with 

25 

respect  to  time 


2 . 1 


The  instantaneous  acceleration  of  a  particular  object  is 

defined  as  the  force  required  to  change  its  velocity  with 

2  5 

respect  to  time 


F 


M 


dV 

dt 
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2 . 2 


When  dealing  with  rocket  systems,  those  forces  must  be 
equal.  When  doing  this,  one  must  note  that  the  change  in 
mass  of  the  propellent  is  the  negative  change  in  mass  of 
the  rocket. 


dm  =  -  dM  2  .  J 


Additionally,  the  propeliant  velocity  can  be  rewritten 
1  6 

as  : 


V 

p 


i 

*  p 


2  •  4 


where  I 

«P 

resul ting 


is  the  specific  impulse  in  units  of  time, 
differential  equation  is: 


The 


dm  _  dV 

M  g  I 


With  M  defined  as  the  payload,  propulsion  system,  and 
power  and  MQ  as  plus  propellant,  tne  resulting  solution 
to  the  differential  equation  is: 


exp 


i 

-  tV  j 


2.6 


The  specific  impulse  is  an  important  element  required  for 
the  velocity  of  the  rocket,  as  seen  from  equation  2.6. 

3  i 

This  specific  impulse  can  be  expressed  in  several  ways 


I 

9  p 


V  P 


f 

fli  g 


g 
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For  a  final  temperature  of  zero  and  isentropic 

expansion,  the  specific  impulse  can  be  calculated  for  pure 

J  .  1  8 

monatomic  or  diatomic  gas 


9  p 


20.9  / 


To 

MW 


( monatomic  ) 


2 . 6 


r 


I 


24.9  / 


To 


(diatomic) 


*  P  MW 

Another  important  aspect  of  specific  impuise  is  for 

calculating  the  power  requirements  of  the  rocket  thruster. 

Using  the  thrust  term  for  force,  the  power  needed  can  be 

calculated  from  the  following  equation: 


2.9 


1 


rok 


F  I 


9  P 


2 . 10 


To  minimize  power  requirements,  the  appropriate 
propulsion  should  be  used  for  its  corresponding  mission. 
For  missions  such  as  planetary  lift-off,  high  thrust  is 
required.  For  long  duration  orbital  or  interplanetary 
travel,  high  specific  impulse  is  required.  As  mentioned 
before,  chemical  propulsion  provides  high  thurst  while 
electric  propulsion  provides  high  specific  impulse. 


2.2.2  Orbit. 


Most  applications  for  electric  rocket  systems  involve 
platform  station  keeping  as  an  orbit  about  a  planet.  Each 
orbit  requires  a  different  velocity.  That  velocity  (  Va ) 
should  be  tangential  to  the  orbit  curvature  and  should  have 
a  magnitude  close  to  the  velocity  resulting  from  the 
gravitational  force  of  the  planet  (  Vg  )  .  In  other  words  for 
the  specific  case  of  a  circular  orbit  pattern,  the 
centrifugal  force  must  equal  gravitational  acceleration. 

As  seen  in  Figure  2.1,  the  object  could  orbit  in  either  an 
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elliptical  or  circular  pattern.  Using  the  nomenclature  of 
that  figure,  the  velocity  equation  needed  by  the  rocket  is: 


_  fZz  f  2  r  r 

1  r  J  e  8  [  r  "  aj  2.11 

The  circular  orbit  is  a  special  case  when  a  =  r. 

2.3  Plasma  Physics. 

The  plasma  state  is  considered  the  "fourth"  state  of 
matter.  As  an  anisotropic  media  in  an  electromagnetic 
environment,  it  is  a  complex  region.  Many  types  of  species 
and  reactions  exist  within  the  plasma  region.  An  attempt 
is  made  to  characterize  those  reactions  u3ing  elementary 
plasma  physics  equations.  Additionally,  an  elementary 
description  of  electromagnetics  is  provided. 

figure  1.1  depicts  some  of  the  discharge  actions  within 
a  piasma,  This  figure  schematically  shows  the  plasma 
region  within  a  microwave  cavity.  The  center  of  the  region 
(I)  is  highly  ionized.  The  next  region  (II)  is  less 
ionized.  The  outer  region  (III)  is  the  recombination 
region  for  excited  molecules.  Excited  species  (FES)  flow 
outward,  while  neutral  species  (FNS)  flow  inward.  The  ions 
are  initially  propagated  by  incoming  microwave  power  (MW). 
Finally,  some  energy  is  dissipated  along  the  walls  of  the 
cylinder  through  heat  by  conduction  /  convection  ( HCC  )  and 
through  radiation  (Rad). 

2.3.1  Reactions . 

Reactions  in  this  media  can  be  viewed  differently  than 
a  typical  chemical  reaction.  For  these  reactions,  the 
propagating  mechanism  is  collision,  with  the  reaction  rate 
being  a  collision  rate.  The  collisions  can  be  divided  into 
four  general  categories  ( eee  Figure  2.2).  These  four 


Neutral-Neutral 


© 

- >  0“) 

© 

Ion-Neutral 

©“> 

^G)~>© 

Radiation 


Figure  2.2  Collision  Types 
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collision  types  contain  the  following  basic  reaction 
26 

mechanisms  .  Disassoc iation  and  recombination  reactions 
for  polyatomic  molecules  are  specific  cases  for  excitation 
and  de -exc i ta t i on  reactions. 


Electron  -  Neutral/Ion 


A  +  e  - 

*  - 

- — *  A  +  e 

electron  excitation 

2.12 

* 

A  +  e - 

- ►  A  +  e 

electron  de-excitation 

2  .  1  J 

A  +  e  - 

- »  A  *  +  2e~ 

electron  ionization 

2.14 

♦ 

A  +  e  - 

- *  A  +  ho 

radiative  recombination 

2.  15 

Ion  - 

Neutral 

A*  +  B  - 

- ►  A  +  B  + 

electron  transfer 

2  .  16 

A+  +  A  - 

- »  2A  +  +  e 

ion  ionization 

2.  17 

Neutral 

-  Neutral 

A  +  A  — — 

* 

- *  A  +  A 

neutral  excitation 

2  .  18 

A  +  A  ' 

- *  A +  +  A  + 

e  neutral  ionization 

2.19 

« 

A  +  B  - 

* 

- *  A  +  B 

excitation  transfer 

2 . 20 

Radiation 

A  +  - 

♦ 

- *  A 

radiative  excitation 

2 .21 

* 

A  - 

- ►  A  +  ho 

radiative  de-exc i tat  ion 

2 . 22 

A*  +  - 

- ♦  A +  +  e 

radiative  ionization 

2 .23 
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2.3.2  Electromagnetic  Field. 


As  implied  by  its  name,  an  electromagnetic  field  is 
composed  of  an  electric  field  (E)  and  a  magnetic  field  (H). 
These  vector  fields  satisfy  four  classic  equations,  known 

30 

as  Maxwell’s  equations  (listed  below) 


?  •  D  =  P 


2 


O 

Li 


4 


7  x  E 

V  x  H 

V  ■  b 


-  dB 

at 

7  ap  d(  tE  ) 

J  *  at  at 

o 


2  .  26 

2 .27 


The  electric  induction  (D)  is  a  function  of  E,  the 
polarization  density  (P),  and  the  permittivity  of  free 
space,  CQ. 


D  =  e  E  +  PQ  2.28 

The  polarization  density  can  be  dropped  from  this  term  if 
we  know  the  electric  susceptibility  (X  )•  The  dielectric 

ft 

constant  (x  )  is  one  plus  X  •  Therefore,  the  electric 
e  e 

induction  can  be  rewritten  as  follows: 

D  =  €  (  1  +  X  )  E  =  €  x  E  =  eE  2.29 

O  e  U  e 

The  magnetic  induction  (B)  is  a  function  of  H,  the 
magnetization  density  IM),  and  the  permeability  of  free 
space,  PQ. 


=  MQ(H  +  M) 


B 


2  .  30 
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Once  again,  the  magnetization  density  can  be  dropped  from 
this  term  if  we  know  the  magnetization  susceptibility  (X  )• 

in 

The  relative  permeability  (x  )  is  one  plus  X  •  Therefore, 

tn  m 

the  magnetic  induction  can  be  written  as  follows: 

B  =  u(X+l)H  =  pxH  =  /iH  2.31 

0  'm  0  m 

The  current  density  iJ)  is  a  scalar  multiple  of  E  plus  the 
velocity  (u)  of  the  media  crossed  with  B.  That  scalar 
function  is  the  electrical  conductivity  ( o )  . 

J  =  «JE+u\B  =  nu  2.32 

e 

The  total  charge  (q)  of  the  media  is  the  total  volume 
integral  of  the  charge  density  (p).  For  the  plasma  state, 
the  q  should  be  zero  because  the  postive  charges  should 
equal  the  negative  ones. 

2.3.3  Conservation  of  Particles. 


The  conservation  of  mass  must  be  upheld.  Matter  can 

neither  be  created  nor  destroyed.  For  this  case,  the 

conservation  of  particles  can  be  defined  as  the 

accumulation  of  electrons  with  respect  to  time  equal  to  the 

negative  divergence  of  the  current  density  plus  the 

27 

generation  of  electrons 

-  *  •  J  +  r  2.33 

e 


<3n 

_ « 

w~ 


dn 

e 


7  •  ( n  u )  +  r 

e  e 


2 . 34 


The  electron  generation  is  composed  of  ionization  minus 

deionization  processes.  That  net  generation  term  can  be 

28. 

written  as  : 
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r=Ptn-<xnn  2.35 

e  l  e  e  i 

For  low  pressures  (roughly  less  than  100  torr) ,  this 
generation  term  is  commonly  neglected.  The  resulting 
equation  with  all  deionization  occuring  on  the  surface  walls 
i  s : 


On 

dt 


-  ?  •  (  n  o  ) 

e 


2 . 36 


The  electron  density  can  be  expressed  as  a  function  of 
position  and  time.  The  plasmas  in  this  experiment  were 
propogated  in  a  cylindrical  cavity.  Using  separation  of 
variables  (radial,  axial,  and  time),  the  electron  density 
can  be  expressed  as2?: 


n  =  R(r)  Z(z)  T(t)  2.37 

e 

For  this  calculation,  the  electron  density  is  assumed  to  be 
at  steady  state.  The  T(t)  term  drops  out.  Assuming  only 
diffusional  decay  and  no  E,  the  neU  term  was  expressed  as 
the  negative  scalar  diffusion  coefficient  (D)  times  the 
gradient  of  the  electron  density.  Then,  equation  2.36 

,  27 

reduces  to  : 


-  V  •  in  u)  *  -  V  •  (-D?n  )  =  DV2n  =0  2.38 

•  e  e 


In  cylindrical  coordinates  (see  Figure  2.3),  this  equation 
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l  s  : 


2  .  .  2  2 

On  On  On  On 

- _  +  - —  +  — - —  +  - =  0  2.3 

dr  r  dr  r  00  0z 

Assuming  radial  symmetry,  the  electron  density  is  not 
dependent  upon  angle.  Thus,  the  partial  of  n^  with  respect 
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to  the  angle  is  zero.  Using  the  separation  of  variables, 

2  ? 

equation  2.39  becomes  : 


d  2  R 

dr2 


dR 

r  dr 


d2Z 


Z  dz 


0 


2 . 40 


A  solution  to  this  equation  can  be  found  by  rewriting 
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equation  2.40  as  : 


j92  +  a2  =  0 


41 


f  2  1 

d  R  dR 

- - -  +  - 

k  dr  r  dr  _ 


2 

-  a 


2 .42 


d  — -  =  -  2.43 

Z  dz 

Equation  2.42  can  be  solved  by  using  substitution  of 
variables27 . 


r 


n 

L* 
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Substituting  equation  2.44  into  2.42  results  in  a  zero 

27 

order  Bessel  equation 


.2  d 2 R 


+  U 


dU 


dR 

dU 


+  U  R  = 


i - 
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Using  the  boundary  condition  that  the  electron  density  is 
zero  at  r  equal  to  the  radius  of  the  cavity,  the  «  terms 

n 

are  the  zeroes  to  the  zero  order  Bessel  function  ( JQ )  of 
the  first  kind.  Thus,  one  can  write  the  solution  to  this 

27 

equation  as  : 
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,v  » 

R(r)  =  Y  B  J  (U  )  =  y  B  J  r) 

Z,  n  0  n  £,  n  O  n 


2 .46 


n  «  1 


n  ■  1 
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The  solution  to  equation  2.43  is  easier  to  find“  . 


Z  (  z  )  =  7  C  Cos ( 0  z ) 

L.  n  n 


2.47 


n  »  1 


Using  the  boundary  condition  that  the  electron  density  is 

.  27 

zero  at  z  =  *1 


Therefore,  the  final  solution  to  this  symmetric  steady 
state  piasma  at  low  pressures  is  trivial.  The  result  is  a 
function  of  a  zero  order  Bessel  function  for  radial 

27 

position  and  a  cosine  function  for  axial  position 


=  n.o2_B»Cn  C°5(  n  2  2  j 


The  constants  B  C  can  be  approximated  to  be  normalized  to 

n  n 

that  of  the  electron  density  at  the  center  (point  0),  with 
the  n@o  being  defined  as  that  electron  density.  These 
terms  can  be  found  by  using  the  properties  of 
orthogonality.  Bessel  functions  and  cosine  functions  are 
series  of  orthogonal  terms.  Thus,  these  terms  can  be 
calculated  as : 


B 

n 


a 

n 


2 _ 

J.  (a  ) 

1  n 


2 . 50 
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C 

n 


4 


n  n 


2.51 


Those  terms  can  be  tabulated  and  are  provided  in  Table  2.1. 


Table  2.1  B  and  C  Values  of  Simple  Electron  Density 

n  n 

Func  t ion 


n 

B 

n 

c 

n 

1 

1 . 602 

1.273 

') 

- 1  .  Oo5 

0 . 000 

3 

0 .851 

-0 .424 

i 

-T 

-0.730 

0 . 000 

2.3.4  Conservation  of  Momentum. 

As  seen  from  tne  simple  model  provided  by  the 
conservation  of  particles,  an  accurate  model  of  the  plasma 
region  can  not  be  accomplished  without  accounting  for 
electron  generation.  This  generation  involves  collision 
physics.  Equations  involving  conservation  of  morentum  are 
needed  to  model  the  generation  term  in  equations  2.34  and 
2.35. 


Classical  scattering  can  be  used  to  describe  the 
momentum  loss  per  collision,  along  with  the  average  number 
of  collisions  expected  within  a  certain  region  (see  Figure 
2.4).  The  average  momentum  loss  per  collision  can  be 
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expressed  by  the  following  : 

M,  =  p  ^  (  1  -  cos  0 )  2.52 

L  me 

The  average  number  of  collision  expected  between  the 
deflection  angle  of  collision  (0)  and  0  +  d0  is2?: 


n  o 

tm 


u;,  s  N  V  2  *  sinS  dfl  2.53 

r  8 

The  average  momentum  loss  can  be  found  by  integrating  the 
momentum  loss  times  the  average  collision  over  all  of  8. 

To  simplify  the  expression,  a  momentum  cross  section  term 

(Q  )  is  defined  as28: 

1 


n 


Q  =  2  n  i  1 

J 

0 


cos  & )  sin  &  d 


2 . 54 


Using  this  cross  section  term,  the  average  momentum  loss 
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can  be  expressed  by  the  following  equation  : 


v  M  >  -  p  v  N  V  $  2.55 

m  ©  r  m 

The  scattering  differential  cross  section  can  be  expressed 

considering  only  coulomb  collisions.  For  this  type  of 

collision,  the  cross  section  can  be  calculated  from  the 

2  8 

differential  element  : 


da 

dO 


e:  e2 


,  2  2 

8  n  c  m  V  sin 
o  e 


ID 


2.56 
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The  collision  frequency  <Uj)  can  now  be  defined^  . 


V.  = 


<.M4  > 


I  u  i> 


N:  to  V  .> 
e 


2.57 


The  excitation  rate  <-<?V  2  can  be  broken  down  into 

e 

specific  excitation  or  de-excitation  terms.  The  following 

z  8 

equation  shouid  hold  true  : 


:0  V  2  =  - 

1 , 0  e  n . 


■<,CT  V  ;> 
0,1  e 


2.58 
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Experimentally,  the  cross  section  term  can  be  found  from 
ion  beam  scattering  using  the  following  expression27; 

Q  *  2  (a,  -  a  In  £  =  2  Q  2.59 

m  1  2  r»  ex 

The  $  is  the  total  cross  section  for  charge  exchange. 

ex 

The  a  term  is  calculated  from  the  linearity  of  Q  versus 
2  ex 

the  log  of  e  where  e  is  the  laboratory  collision  energy. 

n  n 

A  complete  explanation  of  the  theory  and  experimental 
procedure  of  ion  beam  scattering  may  be  found  in  Mason  and 
McDonald’s  book27. 


Deionization  is  needed  to  complete  the  conservation  of 

momentum.  Electron  recombination  is  the  process  of 

deionization.  This  recombination  is  highly  dependant  upon 

the  electron  temperature  (T  ).  For  simplicity,  only 

0 

approximations  will  be  given  for  three  primary  types  of 

27 

recombinations  (radiative,  dissociative,  and  collisional) 
These  values  are  approximated  from  experimental  data  of 
gases  in  various  environments. 
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2.3.5  Conservation  of  Energy. 


Another  important  aspect  of  modeling  the  plasma  region 
is  developing  an  accurate  model  showing  the  temperature 
gradient  within  the  plasma.  The  conservation  of  energy 
equation  can  be  used  for  this  model.  This  energy  equation 
can  be  written  on  a  differential  element  of  the  plasma  at 
the  microscopic  level.  The  accumulation  of  energy  (dE/dt) 
is  equal  to  the  flow  of  energy  in  (Fi )  minus  the  flow  of 
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energy  out.  i  F  )  plus  the  generation  ot"  energy  (r  )  within 

o  en 

9 

the  differential  element  . 


F  -  F  +  r 

i  0  en 


2 .63 


This  energy  can  be  expressed  in  terms  of  power  at  the 
microscopic  level.  The  flow  terms  and  the  generation  term 


can  be  sub-divided  by  the  following  power  terms  : 


F  -Fn=P.+P.+P 

i  0  rad  end  env 
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r  =P+P  +P  -P  -P  -  P 

en  e  roc  eel  el  oxc  ion 


2.65 


The  radiation,  conduction,  and  convection  terras  are 

functions  of  temperature,  heat  capacity,  and  boundary- 

conditions  (such  as  cavity  wall  temperature  and  material 

type)  .  The  electron  recombination  term  may'  be  neglected 

for  low  pressures.  The  electric  field  power  can  be  written 
28 


P  =  a  E 

e 


2.66 


with  the  conductivity  coefficient  written  as' 


n  e 

e 


m  (  5  +  ) 


2.67 


The  excitation  and  ionization  terms  are  functions  of  the 
density  of  the  species  and  of  their  excitation  /  ionization 
rates.  These  two  terms  can  be  written  for  the  electron 


n  n  £  v 

e  n  x  x  e 


2 . 68 


P  =nn£<.ou;  +  nn£  <  o  v  > 

Ion  e  n  i  le  exxi  xie 
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Finally,  the  collision  terms  provide  both  generation  and 

degeneration  of  energy.  The  generation  type  involves 

superelastic  (or  de-excitation)  collisions  and  can  be 
2  8 

express  as  : 


P  ,  =  n  n  £  <*<*£>  > 

se  1  e  x  x  d  © 


2  .  70 


The  degeneration  term  involves  elastic  collisions  of 


28 


electrons  and  neutrai  species  which  can  be  expressed  as  : 


P  =  n  n  (-—-I 

el  e  n  l  M  J 


<a  V  > 

s  e 
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2.4  Discharge  Properties. 


2.4.1  Energy  Distribution 

The  measurement  of  the  energy  distribution  for  this 
research  involved  a  macroscopic  frame  of  reference.  The 
macroscopic  energy  balance  only  covers  energy  entering  and 
leaving  the  entire  system.  The  system  i3  defined  as  the 
microwave  resonance  cavity  (see  Figure  2.4).  An  energy 

9 

balance  written  around  this  system  can  be  written  as  : 

E  =  E+E+E+E  2.72 

a  g  a  w  r 

E#  is  the  energy  entering  the  system  from  the  microwave 
power  source  and  is  represented  by  the  following 
expression : 

E  -  P  -  P  2.73 

a  i  r 

is  tne  energy  absorbed  by  the  air  cooling  and  is 
measured  by  the  air  flow  rate  (F  )  and  its  temperature  rise 
(AT  )  using  the  following  equation  (assuming  the  heat 

a 

9 

capacity,  C  ,  is  constant  over  the  temperature  range)  : 
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E  =  C  F  AT  2.74 

a  p ,  a  a  a 

is  the  energy  absorbed  by  the  water  cooling  and  is 
measured  by  the  water  flow  rate  (F  )  and  its  temperature 

W 

rise  (AT^)  using  an  equation  similar  to  2.74.  Er  is  the 
energy  escaping  the  system  as  emission  radiation.  Because 
the  system  is  almost  entirely  enclosed,  this  radiation  is 
neglected.  Rearranging  the  above  equations  yields  the 
following: 


E  “  P  -  P  -  C  F  AT  -  C  FAT  2.75 

g  1  r  p , a  a  a  p,ww  w 


For  this  research,  the  above  energy  expressions  are 
reported  in  percentages,  as  opposed  to  absolute  values. 
The  following  equations  indicate  how  these  values  are 
computed : 


E 

%E  =  x  100% 

a  E 

8 

E 

%E  =  —=r~  x  100% 

w  E 


2 .76 


2.77 


E 

%E  =  —=r—  x  100% 

a  E 


2.78 


As  determined  from  reading  (measurement)  errors  and 
statistical  deviations,  the  following  absolute  error  range 
estimates  are  provided:  2%  for  both  air  and  water  energy, 
3%  for  power  source  energy,  and  a  composite  of  7%  for  gas 
energy.  These  percentages  are  given  with  respect  to  the 
microwave  source  energy. 


2.4.2  Electromagnetic  Resonator  Modes  in  Microwave  Cavity 

The  experiments  for  this  thesis  use  a  cylindrical 

cavity  resonator.  A  transverse  magnetic  ( TM  )  mode  is 

mnp 
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used.  To  understand  the  electromagnetic  field  in  the 
cavity,  one  roust  analyze  the  three  cylindrical  components 
of  the  electric  (E)  and  magnetic  (H)  fields.  For  a 
circular  TM  mode,  the  axial  (z)  component  of  the 

mnp 

magnetic  field  (H^)  is  zero. 

A  solution  to  Maxwell’s  equations  will  be  calculated 
using  phasors  to  characterize  the  time  aependance  of  the 
wave  fields.  Before  solving  Maxwell’s  equations  for  the 
waveform,  it  may  be  clearer  to  define  a  few  constants 
i terms ) . 


e 


=  -g1  =  ju>  (phasor  constant)  2.80 


k  =  /  k2  ♦  k2 


2.81 


0 


p 


2 .82 


Using  Maxwell’s  equations  (number  2.27  and  2.28)  and 
equation  2.31,  one  could  calculate  the  electric  and 
magnetic  fields  in  the  cavity.  First,  one  must  calculate 
the  waveguide  equations  before  considering  the  resonator. 
Using  the  boundary  condition  that  is  zero,  the  radial 
(r)  and  angular  (9)  components  of  the  electric  and  magnetic 
fields  are  calculated  as  follows  (in  phasor  and  cylindrical 
coord i nates ) 3  * : 


aEz 

r  ae  *  kz  Ee 


j  W  M  H 

r 


k  E 

z  r 


dE 

dr 


j  «  M  H0 


2.83 


2.84 
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k 

Z 

Ho 

=  j  u  f 

E 

r 

2 .85 

k 

2 

H 

r 

=  j  u  e 

Ee 

2 . 86 

These  four  equations  can  be  simplified  by  separating  the 
unknown  field  components  and  rewriting  them  as  functions  of 
one  variable  (E  ).  The  following  are  the  simplified 

33 

equations  : 
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dE 

Z 

k2 

c 

dr 

k 

z 
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2 .87 
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2.89 


2 . 90 


The  solution  to  the  above  equations  could  easily  be 

calculated  if  we  knew  E  .  Using  the  Laplacian  operation 
2 

(V  )  and  separation  of  variables  (radial  and  axial 
functions)  on  E^,  we  can  derive  a  general  solution.  This 
derivation  is  similar  to  that  done  in  section  2.3.3.  The 
following  is  the  general  solution  to  7Ze^  =  0  33 . 


=  C 

m 

L. 


J  (  k  r  ) 

■  C 


BN  !  k  r ) 

mm  c 


ir 

JL 


A  s  i  n  (  )  + 


B  cos ( raO )  exp (  k  z ) 

m  z 


For  finite  values  at  the  origin  ( r  =  0  )  ,  one  obtains  a 

3  3 

reduced  solution  for  E 


E  =  BC  J(kr)cos(ra0)expikz)  2.92 

z  m  m  m  c  z 

E  =  E  J  (k  r)  cos(m6)  exp(-j  P  z) 

Z  m  m  c  p 


.91 


2 . 93 
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Therefore,  the  remaining  five  components  of  this  waveform 
are  : 


^pR0E>T  fP  .nr]  f-n  fPBnr|' 

p  2  rJ4  Ro  J  "  R0  Ro  J 


e  o  s  i  m0  )  expijp  z  i  2.94 

p 


-j  Vo  "  E, 


P  Rn 

m  n  0 


s in ! m© )  exp ( J0p2 ) 
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/P  T'v 
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■jU  CR0E1  "  (P*nr)  jijn 

‘  rJ4  R0  ‘  R0 


6  '  P  2 

m  n 


P  r^, 

cost  m9)  exp(j/3  z)  2.97 
Ro  J  p 


H  =  0 
2 


2.98 


A  table  showing  the  m  zero  of  J  (P  )  follows.  The 

m  mn 

waveform  is  characterized  by  this  zero  and  is  written  as 
TM  . 


Table  2.2  m  Zero  of  J  (P  ) 

n  mn 


n 

m  I 

0  1 

1 

1 

2  1 

3 

11.792 


3 .832 
7.016 
10.173 
13 . 323 


5 

1  36 

8 

417 

1 1 

620 

14 

796 

6 . 380 
9.761 
13.015 
16 . 200 
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For  a  resonator,  one  must  consider  boundary  conditions 

for  E  at  the  base  and  height  of  the  cavity.  Fox  z.-0  and 
z 

z«Lg  (Lg  equals  the  height  of  the  cavity  -  otherwise  known 
as  the  microwave  short  length),  the  following  is  true: 


For  the  radius  at  the  edge  of  the  cavity  ( r«R  )  ,  the 
electric  field  reaches  a  maximum  value. 


bE 


dr 


The  exponential  term  of  Ez  (see  equation  2.93)  may  be 
transformed  to  a  sine  function  of  z  if  we  take  the 
following  value  for  the  propagation  constant: 


The  E  terra  then  reduces  to: 
2 


E  J  ( k  r )  cos ( n@ )  sin 

m  m  c 


p  n  z 


Thus,  the  solution  to  this  equation  defines  a  TM  mode. 

mn  p 

Additionally,  the  frequency  of  the  microwave  must  be  known 
to  calculate  k  (see  equations  2.79  through  2.82)  .  The 

C 

frequency  (u)  is  2.45  GHz.  This  frequency  also 
characterizes  the  electric  and  magnetic  field  components. 
For  this  experiment,  TM  1  and  TM012  modes  were  used  (see 
Figure  2.5  for  wave  patterns  of  those  modes). 


2.4.3  Spectroscopic  Analysis. 

39 

Free  electrons  are  an  integral  part  of  a  plasma 
Thus,  ionization  and  deionization  processes  within  a  plasma 
occur  frequently  and  may  be  thought  of  as  the  plasma's 
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'igure  2.6  Experimental  TK  Modes 
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sustainment  (or  driving  force).  As  such,  the  electron 
temperature  can  be  seen  as  one  of  the  few  properties 
characterizing  the  state  of  a  plasma.  One  way  to  measure 
the  electron  temperature  is  through  spectroscopy. 

There  are  two  distinct  methods  for  conducting 
spectroscopic  measurements,  absorption  and  emission.  For 
this  research,  emission  spectroscopy  is  used.  The 
spectrometer  detects  the  spectral  lines  emitted  by  the  gas 
as  the  atom  decays  from  the  excited  state  to  the  ground 
state.  In  the  plasma,  the  atoms  are  excited  from 

38 

electromagnetic  radiation  anu  from  collisions  of  species 


The  intensity  of  the  spectral  line  is  determined  by 
two  factors,  the  Boltzmann  distribution  and  the 
transmission  moment  (with  Einstein  coefficients).  The 
Boltzmann  distribution  law  characterizes  the  population  of 
various  excited  levels  (N  )  using  the  following  equation 

”  38 

i with  the  degeneracy  terra,  g  ,  added)  : 

n 


N 

n 


g  N 


2  .  103 


The  transition  moments  are  determined  through  the 
interactions  of  the  species  with  the  emitted  radiation. 

Using  the  relative  emission  line  intensities  for 
spectroscopic  measurements,  an  emission  line  intensity  going 

g 

through  the  spectrometer  can  be  given  by  : 


I  = 
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N  h  g  A  R.  v 

n  nm  A.  n» 


dQ 


4n 


2.104 


For  ground  state  transitions  when  Q  becomes  independent  of 
temperatures,  this  equation  can  be  reduced  to  provide 
an  easier  equation  to  calculate  the  electronic  temperature. 
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For  this  thesis,  local  thermal  equilibrium  is  assumed. 
Although  this  assumption  may  not  provide  accurate  resuits, 
it  does  allow  for  a  quick  "ball- park"  calculation  of  the 
electron  temperature.  This  assumption  allows  us  to  assume 
that  the  electron  temperature  is  approximately  equal  to  the 
electronic  (excitation)  temperature. 


T  T  2.106 

e  elc 

Errors  are  involved  in  this  techinique  to  measure  the 

electronic  temperature.  Although  we  can  easily  measure  the 

relative  population  of  an  excited  species,  the  absolute 

value  is  not  accurately  obtained.  This  comes  from  the 

procedure  of  measuring  relative  einision  line.  An  absolute 

emision  line  can  not  be  measured  in  this  spectroscopy 

experiment.  Additionally,  emission  spectroscopy  does  not 

provide  information  about  the  ground  state  species,  which 

happens  to  be  the  most  populated  state.  According  to  an 

analysis  by  Chapman,  we  can  expect  errors  in  excess  of  50% 

0 

because  of  these  limitations  . 


A  better  way  of  measuring  the  electron  temperature  is 
measuring  line  width  species  such  as  hydrogen  which 
undergoes  a  linear  Stark  broadening  and  which  does  not 
assume  local  thermal  equilibrium. 


CHAPTER  III 


ELECTRIC  PROPULSION 


3.1  Characteristics. 

Electrical  propulsion  systems  use  charged  particles 
accelerated  by  an  electric  field  as  a  working  fluid.  These 
systems  are  capable  of  creating  greater  specific  impulses 
than  chemical  or  nuclear  propulsion  systems.  Likewise,  this 
capability  further  supports  their  suitability  for  steering 
rockets  in  space.  There  are  three  basic  types  of  electric 
rocket  thrusters:  electrothermal,  electrostatic,  aid 
electromagnetic* . 

Electrothermal  thrusters  use  electric  energy  to  power 
an  arc  or  resistance  heater  to  heat  a  conventional  working 
fluid.  Ions  or  colloidal  particles  make  up  the  working 
fluid  in  an  electrostatic  system  and  are  accelerated  by  an 
electrostatic  field.  Finally,  a  travelling  magnetic  and 
electric  field  system  accelerates  a  plasma  in  an 
electromagnetic  system. 

3.2  Propellants. 

As  mentioned  in  chapter  two,  specific  impulse  is  an 
important  parameter  for  propulsion  systems.  A  higher 
specific  impulse  means  lower  propellant  flow  rate  to  produce 
a  given  thrust  (see  equations  2.1  and  2.4).  In  comparison, 
chemical  systems  have  specific  impulses  in  the  range  of  250 

-  450  seconds,  while  electric  systems  have  the  range  of  300 

-  5000*  seconds20. 
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On  the  negative  side,  electric  propulsion  systems  are 
characterized  by  low  thrust20.  In  other  words,  these 
systems  do  not  possess  the  force  needed  to  quickly  overcome 
a  strong  gravitational  counteraction.  This  requires  longer 
operations  to  achieve  the  desired  velocity  change  than 
chemical  propulsion.  For  example,  a  mission  may  require 
several  hours  of  operation  per  day  for  electric  propulsion. 
An  equivalent  mission  for  chemical  propulsion  would  require 
several  minutes  per  week. 

Why  is  propellant  amount  so  i rapo rtant ?  The  propellant 
is  a  major  factor  in  the  cost  of  a  spacecraft  mission.  In 
one  example,  the  propellant  accounted  for  43%  of  the  mass 
in  the  Galileo  mission.  In  another  example,  76%  of  the 
mass  in  the  Comet  Rendezvous  Asteroid  Flyby  (CRAF)  mission 
was  propellant. 

3.3  Applications. 

Currently,  electrothermal  propulsion  performs 

Nortn-South  Stat ionKeeping  (NSSK)  on  many  geosynchronous 
16  2  1 

satellites  '  .  American  Telephone  &  Telegraph  Company 

(AT&T)  include  electric  propulsion  options  on  their  newest 
satellites  -  Telstar  IV.  This  type  of  propulsion  will  be 
used  on  the  Space  Station  Freedom  as  well  as  other  man 
occupied  platforms16'21. 

Figure  3.1  schematically  illustrates  a  basic  NSSK 
satellite.  The  solar  panels  (power  source)  define  the 
North-South  axis.  For  a  NSSK,  the  rocket  thrusters  act  in 
the  direction  of  the  solar  panels.  This  thruster 
configuration  is  more  efficient  than  East-West 
Stat ionKeeping  (EWSK).  Typically,  two  pairs  of  rocket 
thrusters  are  used.  These  pairs  are  located  at  the  north 
and  soutn  side  of  the  satellite.  Although  only  one 


Figure  3«-  NSSK  Satellite 
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thruster  is  needed  for  each  side,  two  are  used  with  one  as 
a  backup  in  case  the  other  malfunctions. 

Propellant  contamination  is  also  a  concern  for 
application  purposes.  Problems  arise  when  communication 
must  be  transmitted  through  some  portion  of  the  thruster 
plume.  The  radio  frequency  (RF)  signal  may  interact  with 
charged  species  of  that  plume.  This  interaction  is  serious 
as  the  following  r.iay  occur:  reflection  of  the  signal, 
attenuation  and  chase  shift  of  the  signal,  or  generation  of 
noise  in  the  signal. 

3.4  Proposed  System. 

One  type  of  e J ectrotherma]  propulsion  system  uses  a 

6 

microwave  induced  plasma  .  Although  this  system  uses  an 
electromagnetic  environment,  it  is  classified  as  an 
electrothermal  because  it  uses  a  nozzle  (not  the 
electromagnetic  field)  to  accelerate  the  propellant. 
Schematically  shown  in  Figure  3.2  is  a  version  of  this 
system . 

The  power  is  beamed  to  the  spacecraft  from  an  outside 
source  (such  as  a  space  station  or  planetary  base)  as 
microwave  or  millimeter  wave  power.  This  power  is  focused 
onto  a  resonant  cavity  to  sustain  a  plasma  in  the  working 
fluid  via  heating  it.  The  hot  gas  would  expand  through  a 
nozzle  to  produce  thrust.  Alternately  in  a  self-contained 
situation,  power  from  solar  panels  or  from  nuclear  reaction 
could  be  used  to  run  a  microwave  frequency  oscillator  to 
sustain  the  plasma. 
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CHAPTER  IV 


EXPERIMENTAL  SYSTEM 


4 . 1  Introduction. 

The  system  used  in  this  experiment  was  designed  to 
conduct  diagnostic  measurements  of  three  elements  of  plasma 
characteristics.  At  the  macroscopic  level,  the  power 
distribution  and  plasma  dimensions  were  determined  using 
thermocouples  and  visual  photography  respectively.  And  at 
the  microscopic  level,  the  electron  temperature  were 
measured  using  an  optical  emission  spectrometer.  See  Figure 

4.1  for  the  overall  set-up. 

4.2  Microwave  Cavity. 

An  electromagnetic  system  was  needed  to  generate  a 
plasma.  The  microwave  cavity  body  was  made  from  a  17.8  cm 
inner  diameter  brass  tube.  As  seen  in  Figure  4.2,  the 
cavity  contained  a  sliding  short  and  a  coupling  probe  (the 
two  major  mechanical  moving  parts  of  the  cavity).  The 
movement  of  this  short  allowed  the  cavity  to  have  a  length 
varying  from  6  to  16  cm.  The  coupling  probe  acted  as  an 
antenna  which  transmitted  the  microwave  power  to  the  cavity. 
The  sliding  short  and  coupling  probe  were  adjusted  (or 
moved)  to  obtain  the  desired  resonant  mode.  A  resonant  mode 
represents  an  eigenvalue  of  the  solution  to  Maxwell’s 
equations.  Two  separate  resonance  modes  were  used  in  these 
experiments:  TM  (L  -  7.2  cm)  and  TM„,_  (L  =  14.4  cm)  . 

o 1 1  i  012* 

Additional  features  of  this  cavity  included:  two  copper 
screen  windows  located  at  90  degree  angles  from  the  coupling 
probe  (which  allowed  photographic  and  spectral 
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Figure  4.1  Experimental  Set-up 
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measurements),  and  two  circular  holes  (in  bot;i  the  base  and 
top  plates)  to  allow  propellant  and  cooling  air  flows 
through  the  cavity. 

4.3  Plasma  Containment. 

The  plasma  was  generated  in  quartz  tubes  placed  within 
the  cavity  (see  Figure  4.3).  The  inner  tube  is  33  mm  outer 
diameter  and  was  used  for  the  propellant  flow.  The  outer 
tube  was  50  mm  outer  diameter  and  was  used  for  air  cooling 
of  the  inner  tube.  Both  tubes  were  about  2  1/2  feet  long 
and  was  epoxied  to  aluminum  collars.  These  collars  fed  the 
gas  and  air  to  and  from  the  cavity.  For  additional 
protection,  water  cooling  was  done  on  the  collar  downstream 
of  the  cavity. 

4.4  Flow  System. 

Flow  of  99.99%  pure  nitrogen  and  helium  was  controlled 
using  a  back  pressure  regulator  and  a  3/4  inch  valve  in 
front  of  the  vacuum  pump.  A  Heise  gauge  with  a  range  from 
1-1600  torr  was  used  to  measure  the  pressure  of  the  plasma 
chamber.  Four  sets  of  flow  meters  were  used  to  measure  the 
gas,  water,  and  air  flows.  Thermocouples  were  used  to 
measure  the  temperature  of  the  air  and  water  both  entering 
and  exiting  the  cavity. 

4.5  Microwave  Power. 

A  Micro-Now  420B1  (0-500  watt)  microwave  power 
oscillator  was  used  to  send  up  to  400  watts  of  power  at  a 
fixed  frequency  of  2.45  GHz  to  the  cavity  (see  Figure  4.4). 
Although  rated  for  500  watts,  energy  was  lost  f roc  the 
microwave  cable,  circulator,  and  bidirectional  coaxial 
coupler . 


Gas  Outlet 
Water  Inlet/Outlet 
Air  Outlets 


Figure  4.4  Microwave  Power  Source 


Connected  to  the  microwave  oscillator  was  a  Ferrite 
2620  circulator.  This  circulator  provided  at  least  20  dB  of 
isolation  to  each  the  incident  and  reflected  power  sensors. 
The  circulator  protected  the  magnetron  in  the  oscillator 
from  reflected  signals  and  increased  the  accuracy  of  the 
power  measurements.  The  reflected  power  was  absorbed  by  the 
Termaline  8201  coaxial  resistor.  The  incident  and  reflected 
powers  were  measured  using  Hewlett-Packard  8481a  power 
sensors  and  4J5A  power  meters. 

4.6  Temperature  Probes. 

Type  T  thermocouples  (copper  constantan)  with  braided 
glass  insulation  were  placed  at  the  inlet  and  outlet  for  the 
water  and  air  cooling  (see  Figure  4.1).  An  Omega  400B 
Digicator  was  used  to  measure  the  temperature  at  these  four 
locations. 

4.7  Spectroscopy. 

The  radiation  emitted  by  the  plasma  was  measured  using 
a  McPherson  Model  216.5  Half  Meter  Scanning  Monochromator 
and  photomultiplier  detector.  A  high  voltage  of  900  volts 
was  provided  to  the  photomultiplier  tube  ( PMT )  using  a 
Harrison  (Hewlett-Packard)  Model  6110A  (DC)  power  supply. 

The  output  from  the  PMT  was  processed  through  a  Keithly 
Model  616  digital  electrometer.  The  processed  output  is 
sent  to  a  Metrabyte  data  acquisition  &  control  system  and 
recorded  on  a  Zenith  80286  personal  computer  (see  Figure 
4.5). 


The  monochromator  was  positioned  about  100  cm  from  the 
plasma.  The  emission  radiation  was  focussed  on  the 
monochromator  using  two  25  cm  focal  length  glass  lenses. 

This  lens  system  concentrated  the  emission  radiation  on  the 
entrance  slit  opening  of  the  monochromator.  To  optimize  the 


figure  4.5  Spectroscopy  System 


4  B 


intensity  of  the  spectroscopic  emissions,  the  slit  widths 
for  this  experiment  were  set  at  100  microns  for  the  entrance 
siit  and  50  microns  for  the  exit  slit.  The  atomic  spectra 
was  taken  using  the  1200  grooves  per  mm  grating  (plate)  with 
a  range  of  1050  -  10000  A.  This  groove  setting  allowed  for 

a  large  range  of  wavelengths  to  be  observed.  The  reciprocal 

■> 

linear  dispersion  was  16.6  A  per  mm.  The  focal  length  of 
the  spectrometer  was  one  half  meter. 


CHAPTER  V 


ENERGY  DISTRIBUTION 


5.1  Pressure  Dependence. 

The  energy  distribution  within  the  microwave  cavity 

was  analyzed  over  various  pressures,  ranging  from  200  to 

800  torr.  Using  the  same  microwave  system,  Hoekstra  had 

conducted  this  experiment  using  the  two  separate  modes 

(TMotl  and  TMQ 1 2 i  and  the  two  pure  gases  (nitrogen  and 

helium)  .  In  each  of  these  experiments,  the  power  used 

was  about  250  watts  with  air  cooling  flow  of  2  SCFM  and 

water  cooling  flow  of  5.75  mi/sec.  Hoekstra’s  experiments 

involved  gas  flow  in  the  direction  of  the  gravitational 

force  (F  ).  An  additional  experiment  was  done  reversing 
9 

the  flow  of  the  gas.  Using  this  new  data  and  Hoekstra’s 
results  for  the  same  experiment,  an  estimate  is  done  to 
determine  the  gravitational  effects  on  the  other  three 
Hoekstra  pressure  experiments. 

Reverse  of  the  gas  flow  showed  a  slight  change  for 

energy  absorbed  by  the  gas  (see  Figure  5.1),  Nevertheless, 

a  significant  difference  in  the  power  distribution  was 

observed.  About  5%  of  the  total  energy,  which  was  absorbed 

by  the  cavity  wall,  was  redistributed  to  the  air  cooling. 

Thorough  cleaning  of  the  cavity  before  the  new  experiment 

may  be  the  main  cause  of  this  difference.  Dirt  and  oil 

deposits  on  the  cavity  walls  would  absorb  a  large  amount  of 

input  energy.  Thus,  removing  this  dirt  and  oil  wouid 

account  for  the  large  energy  distribution  change  from  the 

wall.  This  new  energy  distribution  agrees  with  the  work 
0 

done  by  Chapman  .  Figures  5.2  -  5.4  shows  estimated  energy 
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Power  Distribution  Versus  Pressure 

(Helium  -  TM012  mode) 
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Power  Distribution  Versus  Pressure 

(Nitrogen  —  TM012  mode) 
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distribution  values  for  the  reversed  gas  flows  compared  to 
Hoekstra . 

5.2  Flow  Dependence. 

The  power  absorption  by  the  helium  gas  was  calculated 
over  a  range  of  various  flow  rates  at  a  constant  pressure 
of  200  torr  using  an  air  cooling  rate  of  2  SCFM  and  water 
cooling  rate  of  5.75  ml/sec.  Tins  experiment  was  oniy 
conducted  in  the  TMoi2  mode.  Once  again  using  Hoekstra’s 
results,  an  estimate  is  made  for  each  mode  (see  Figure 
5.5).  Hoekstra  conducted  a  similar  experiment  in  the  TM 
mode  . 


With  the  exception  of  the  power  distribution  (which 
may  be  caused  by  dirt  on  the  cavity  wall),  the 
gravitational  force  had  an  insignificant  effect  on  the 
plasma .  The  difference  between  the  data  for  flows  with  and 
against  gravity  were  within  the  range  of  experimental  error. 
A  difference  caused  by  gravity  may  be  observed  by 
significantly  reducing  this  error.  Unfortunately,  that 
difference  can  not  be  observed  using  the  current 
experimental  system. 

5.3  Plasma  Power  Absorption. 

The  power  initially  absorbed  by  the  plasma  is  the 
combination  of  the  power  absorbed  by  the  gas  and  by  the 
cooling  air.  For  helium,  power  absorbed  by  the  plasma  was 
about  67%  in  the  TMol2  mode  and  about  80%  in  the  TMQ 1 1  mode. 
The  remaining  power  was  that  power  which  is  absorbed  by  the 
cavity  walls  through  radiation  and  convection. 

This  plasma  power  absorption  was  the  initial 
anticipated  power  we  can  expect  to  use  for  propulsion.  As 
an  example,  Figure  5.6  shows  the  total  power  absorbed  by  the 


Gas  Power  Absorbtion  Versus  Flow  Rate 
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Figure  5.6  Calorimetry  Graph  VI 
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plasma  for  various  flow  rates,  which  was  dramatically 
different  than  the  power  absorbed  by  the  gas  as  shown  in 
Figure  5.5.  The  power  absorbed  by  the  wall  can  be  recovered 
by  such  methods  as  using  heat  exchangers.  A  heat  exchanger- 
can  transfer  energy  to  the  propellant  before  it  enters  the 
cavity. 


CHAPTER  VI 


PLASMA  DIMENSIONS 


6.1  Introduction. 

An  important  element  of  microwave  generated  plasmas  is 
its  size  and  shape.  Its  size  varies  about  pressure,  fiow, 
and  power  changes.  These  measurements  can  be  used  to  model 
heat  transfer  and  chemical  processes  within  a  plasma. 

Plasma  dimensions  were  calculated  from  35  mm  color  slide 
pictures.  The  projected  images  were  measured  and  compared 
with  a  known  grid  pattern.  This  grid  pattern  is 
photographed  and  projected  at  the  same  distance  as  the 
piasma.  Water  cooling  flow  rate  of  5.75  ml/sec  and  air 
cooling  flow  rates  of  2  SCFM  for  helium  gas  and  3  SCFM  for 
nitrogen  gas  were  used. 

A  35  mm  camera,  on  a  tripod,  was  position  about  2  cm 
from  the  microwave  cavity  window  {see  Figure  4.2).  Pictures 
were  taken  using  a  50  mm  lens,  using  Ektachrome  200  ISO 
color  slide  film.  Several  aperture  and  shutter  speed 
settings  were  analyzed  to  determine  the  optimal  settings. 

The  following  data  were  reported  for  aperture  and  shutter 
speed  settings  of  f2.8  and  l/250s,  respectively.  Using  only 
the  TMQ1?,  mode,  four  sperate  measurements  were  made  for  each 
plasma:  height  and  width  for  strong  and  weak  regions. 

Plasmas  for  both  helium  and  nitrogen  were  observed  to  have 
two  distinct  regions  -  a  more  intense  lighter  color  inner 
region  surrounded  by  a  diffuse  darker  outer  region.  Thus, 
it  was  appropriate  to  record  two  separate  measurements , 
identifying  the  inner  and  outer  regions  as  strong  and  weak 
ionization  regions  respectively. 
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5  9 


6.2  Pressure  and  Flow  Dependency. 

Measurements  were  done  on  helium  plasmas  in  the 
pressure  range  of  400  -  1000  torr  using  a  power  source  of 
about  250  watts.  Three  gas  flow  rates  (0,  572,  and  1144 
SCCM )  were  used  for  each  pressure  point  (see  Figures  6.1  - 
6.3).  Nitrogen  was  observed  in  the  pressure  region  of  400  - 
500  torr  using  the  same  power.  Three  gas  flow  rates  (0, 

102,  204  SCCM )  were  used  for  each  pressure  observation  (see 
Figures  6.45  -  6.6).  It  was  noticed  that  plasma  Bize 
decreases  with  increases  in  pressure  and  increases  slightly 
with  increases  in  gas  flow  rate. 

6.3  Power  Dependency. 

Measurements  were  done  on  helium  plasmas  at  a  constant 
pressure  of  400  torr  with  no  gas  flow  for  both  nitrogen  and 
helium.  The  power  was  varied  between  100  -  250  watts  for 
helium  (see  Figures  6.7  -  6.8)  and  between  210  -  260  watts 
for  nitrogen  (see  Figure  6.9).  The  primary  observation  was 
that  small  changes  in  power  at  these  pressure  ranges 
affected  the  plasma  dimensions.  Increases  in  power  resulted 
in  a  Larger  plasma,  as  expected. 

6.4  Plasma  Color  and  Shape. 

For  both  helium  and  nitrogen  plasmas,  the  strong 
ionization  regions  displayed  an  intense  white  color.  The 
color  difference  is  observed  in  the  weak  ionization  region. 
Helium  (displayed  a  purple  color,  while  nitrogen  displayed  an 
orange  ccloi.  Additionally,  these  two  gases  displayed  a 
different,  shaped  plasma.  Both  shapes  can  simply  be 
represented  by  an  oblate  ellipsoid.  The  helium  plasma 
differed  slightly  in  the  middle  with  a  small  indentation, 
resulting  in  a  "dumbbell  shaped"  figure.  Figure  6.10  shows 


Helium  Plasma  Dimensions  versus  Pressure 

(0  SCCM  Gas  Flow  -  TM  012  mode) 
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(572  SCCM  Gas  Flow  -  TM  012  mode) 
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P’igure  6.2  Plasma  Dimensions  Graph  II 
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Helium  Plasma  Dimensions  versus  Pressure 

(1144  SCCM  Gas  Flow  -  TM  012  mode) 
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Figure  6.5  Plasma  Dimensions  Graph 


Nitrogen  Plasma  Dimensions  versus  Pressure 

(0  SCCM  Gas  Flow  -  TM  012  mode) 
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Plasma  Dimensions  Graph  IV 


Nitrogen  Plasma  Dimensions  versus  Pressure 

(204  SCCM  Gas  Flow  -  TM  012  mode) 


Helium  Plasma  Dimensions  versus  Power 

(Strong  ionization  -  TM  012  mode) 
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a  photograph  of  a  nitrogen  plasma  in  a  TMqi2  mode  (courtesy 
of  Maris  Mantenieks  of  NASA  -  Lewis  Research  Center). 

6.5  Plasma  Volume. 


The  plasma  volume  was  easily  obtained  assuming  axial 
symmetry  and  neglecting  the  indentation  of  the  dumbbell 
shape.  The  volume  was  calculated  as  though  the  plasma 
region  behaved  as  an  oblate  ellipsoid.  The  simple  equation 
for  this  calculation  is: 


Volume 


4rt 

Height  ^ 

X 

Width  1 

1  X 

l  2  J 

l  2  J 

6  .  1 


The  data  from  Figures  6.1  -  6.6  were  calculated  for  this 
volume  determination.  This  calculated  volume  data  is 
provided  in  Table  6.1. 


As  seen  from  this  data,  the  volume  was  dependant  upon 
both  pressure  and  flow  rate.  When  pressure  decreased  or 
flow  rate  increased,  the  volume  of  the  plasma  increased. 

6.6  Mechanical  Observations. 


A  small  mechanical  measuring  device  was  constructed  by 
Hoekstra  and  modified  by  Haraburda  to  measure  objects  from 
a  distance.  As  shown  in  Figures  6.11  and  6.12,  an  8  x  7 
inch  object  constructed  of  two  steel  plates  with  a  fixed 
hole  on  one  side  and  an  adjustable  iris  on  the  other  was  a 
device  designed  to  provide  immediate  feedback  information 
concerning  plasma  dimensions.  Using  this  device,  one  could 
obtain  a  reasonable  measurement  within  minutes.  On  the 
other  hand,  a  more  accurate  measurement  could  be  obtained 
through  photography  within  hours.  Thus,  this  device  could  be 
used  during  an  experiment  to  determine  any  irregularities 
in  the  dimension.  This  was  helpful  because  the  experimental 
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Table  6.1  Helium  and  Nitrogen  Plasma  Volumes 
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figure  6.11  Mechanical  Measuring  Device  - 
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run  was  normally  terminated  before  the  photographic 
measurement . 

To  use  this  device,  the  iris  side  is  placed  against 
the  microwave  viewing  window.  The  person  looks  through  the 
other  hole  at  the  plasma  and  adjusts  the  iris  to  contain 
the  plasma,  both  its  height  and  width.  Using  the  scale 
provided,  a  plasma  dimension  can  be  obtained  by  multiplying 
the  reading  value  by  a  calibrated  number.  This  device  has 
a  measurement  error  of  10%. 


CHAPTER  VII 


SPECTROSCOPY 


7 . 1  Introduction . 

The  electron  temperature  was  an  important  parameter 
characterizing  the  plasma.  The  spectroscopy  experiment  was 
used  to  calculate  (or  approximate  a  calculation  of)  that 
temperature.  Chapter  two  described  the  theory  behind  this 
calculation  and  chapter  four  described  the  experimental 
system . 

Before  doing  spectroscopic  experiments,  two  things  were 
done.  First.,  a  computer  program  was  written  (in  Bosic 
language)  to  use  the  data  acquisition  system.  The 
adjustment  coefficients  in  the  program  and  the  frequency 
(timing)  of  data  acquisition  were  optimized  such  that  the 
measurement  errors  were  less  than  0.5%  from  the  digital 
electrometer . 

Second,  the  spectrometer  was  calibrated  using  a  hodei 
245C  tungsten  lamp  (from  Optronics  Laboratories,  Inc.  )  with 
known  intensity  readings.  The  known  intensity  readings  were 
obtained  from  the  17  AUG  82  calibration  of  the  lamp.  A 
current  of  6.500  1  0.001  amps  was  supplied  to  the  lamp. 
Instrumental  readings  were  taken  for  the  spectrometer  over  a 

o  c* 

wavelength  range  of  3000  A  through  9000  A.  The  results  of 
this  calibration  are  provided  in  Table  7.1.  Linear 
interpolation  was  used  to  calculate  the  spectral  response 
function  ( )  for  specific  wavelengths.  The  ; s  in  Table 
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7.1  were  calculated  using  the  following  non-dimensional ized 
equations : 


R, 


Measured  Intensity 
^  Known  Intensity 


=  120 


Known  Intensity 
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Known  Intensity 


Measured  Intensity 


Measured  Intensity 
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To  use  equation  2.105,  several  coefficients  must  be 
provided.  Because  the  electron  temperature  was  calculated 
from  the  slope  of  the  equation  calculations,  several 
wavelengths  must  be  measured.  From  several  scans  of  the 
spectrum,  four  strong  and  recogn i zeabl e  transition  regions 
(wavelengths)  were  observed.  The  data  for  those  wavelengths 
are  provided  in  table  7.21.  Although  equation  2.105 
requires  those  transitions  to  be  ground  state  originating, 
not  enough  strong  transition  regions  were  observed  to 
calculate  an  electron  temperature.  Therefore,  equation 
2.105  was  used  for  non-ground  state  transitions. 

Table  7.2  Electronic  Transition  Values  for  Helium 
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E 
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0.087 

9 

0.0320 
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3888.65 
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7.2  Experimental  Results. 

An  experiment  was  done  using  helium  in  the  TM  mode 
with  no  gas  flow  at  a  power  of  220  watts.  The  pressure 
range  was  400  -  800  torr.  For  this  experiment,  the  electron 
temperature  was  assumed  to  be  that  of  the  electronic 
temperature  under  the  assumption  that  local  thermal 
equilibrium  occurs  (see  equation  2.106).  The  results  are 
provided  in  Figure  7.1. 

Because  there  was  a  shape  and  magnitude  difference 
between  my  data  and  that  of  Hoekstra,  I  plotted  the  same 

3  2 

data  for  a  similar  experiment  from  Mueller  and  Micci 
Although  not  shown,  Chapman's  results  were  important  in  that 
his  results  showed  the  electron  temperature  dependence  upon 
pressure  .  Add itionally,  it  shows  results  using  Hoekstra’s 
technique  of  data  acquisition  and  interpretation. 

The  volume  electron  temperature  is  expected  to  decrease 

as  pressure  increases.  It  takes  less  energy  to  maintain  a 

smaller  volume  of  plasma  (which  decreases  with  increasing 

pressure).  However,  the  peak  electron  temperature  is 

expected  to  increase  with  pressure  increases.  Chapman's 

results  for  very  low  pressure  (0.5  -  10  torr)  showed  a 

a 

decreasing  temperature  with  an  increasing  pressure  . 

The  magnitude  difference  can  be  explained  because  a 
different  technique  in  data  acquisition  occurred.  Hoekstra 
used  the  same  technique  Chapman  did  by  measuring  the  peak 
heights  on  a  chart  recorder.  Chapman’s  results  had 
temperatures  around  5,000  K.  My  data  were  calculated 
measuring  tne  area  under  the  peak  using  a  data  acquisition 
system  with  a  computer. 


ELECTRON  TEMPERATURE 
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Therefore,  the  electron  temperature  is  expected  to  be 

o 

near  13,000  K  for  pressures  near  atmospheric.  Additionally, 
the  electron  temperature  is  expected  to  decrease  with 
increasing  pressure. 


CHAPTER  VIII 


COMPUTER  MODEL  OF  ELECTRON  DIFFUSION 


8.1  Introduction. 


An  explanation  of  the  diffusion  of  the  electron  in  the 
plasma  region  was  briefly  covered  in  chapter  two.  An 
analytical  solution  was  obtained  by  neglecting  the 
generation  term  in  the  continuity  equation.  As  previously 
mentioned,  this  is  typically  assumed  for  very  low  pressures. 
Unfortunately,  this  assurapt ion  is  not  valid  for  high 
pressures,  such  as  atmospheric.  As  pressure  increases,  the 
density  of  the  gas  increases.  This  increased  density 
results  in  an  increase  in  the  collision  rates,  such  as 
recombination  and  ionization.  The  purpose  of  this  chapter 
is  to  develop  a  numerical  method  to  see  the  effects  that 
electron  recombination  and  electron  ionization  have  upon  the 
distribution  gradient  of  the  electron  density  within  the 
plasma . 

8.2  Development  of  Mathematical  Model. 


An  assumption  is  made  that  the  ion  and  electron 
densities  are  equal.  As  a  result,  equation  2.34  (continuity 
equation  for  electrons)  can  be  re-written  as  (with  N  defined 
as  the  normalized  electron  density): 


d2  N  +  dN  d2N 

drZ  r  dr  d z2 


a 


N  2  -  i>  (  N 
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This  equation  includes  ambipolar  diffusion.  For  no 
generation,  this  equation  can  be  simplified  by  separating  N 
into  separate  variables  as  mentioned  in  chapter.  However, 
this  method  will  not  work  for  generation  because  of  the 
non-linear  recombination  term  on  the  right  hand  side.  This 
tnesis  oniy  considered  the  radial  part  of  this  equation.  To 
calculate  the  electron  density  along  the  radial  axis  about 
the  center  plane  ( z=0 j  can  be  done  by  saying  that: 


2 

02 


Tiiis  is  valid  because  the  change  of  N  is  symmetric  about  the 
origin.  The  resulting  equation  is  a  second  order 
differential  equation: 


This  equation  was  re-written  as  a  set  of  two  first  order 
differential  equations  by  setting: 


8*4 


Therefore,  the  following  set  of  first,  order  differential 
equations  were  used  to  determine  the  electron  density 
distribution  in  the  radial  direction  along  the  center  plane : 
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These  equations  were  solved  using  the  normalized  density. 

The  normalized  density  was  the  non-d  imeus  i  onai  quantity  with 
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N  equal  to  one  when  the  electron  density  is  at  its  maximum 
value  (or  the  center  of  the  plasma). 


For  this  simulation,  the  maximum  density  wili  be 

This  density  uas  used 


12  3 

assumed  to  be  1x10  electrons  /  cm 


because  Chapman  observed  electron  densities  around  this 
value'1.  Using  this  assumption,  the  recombination 
coefficient  must  be  normalized.  This  coefficient  is  a 
function  of  temperature.  For  the  assumed  maximum  electron 
density  and  for  a  temperature  range  of  250  -  64,000  K,  the 


recombination  coefficient  has  a  range  of  1x10 

,  .  3  4  0 

recombinations  cm  /  sec 


-  1  3 


1x10 


The  two  boundary  conditions  required  to  solve  this 
problem  are : 


at  r  =  R 
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N  = 


at  r  =  0 
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8.3  Numerical  Analysis. 


A  variable-step  Runge  Kutta  method  will  be  used  to 
solve  these  two  differential  equations  with  the  given 
boundary  conditions.  A  fourth  and 

fifth  order  evaluation  with  six  evaluations  per  step  will  be 
done.  These  two  evaluations  will  be  done  to  estimate  the 
relative  error  of  each  step.  The  following  formulas  for 


solving  equation  8.5  will  be  used  in  the  computer  program 
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The  above  formulas  are  the  save  for  solving  equation  8.6. 
Because  those  two  equations  (8.5  and  8.6)  are  coupled,  they 
must  be  solved  simultaneously.  The  relative  error  per  unit 
step  can  be  written  as: 


N (  e  '  1  ’  -  N  |  5  e  INC — - 

itl  * ♦ 1 1  rellb-a) 
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The  variable  step  in  this  program  is  the  recalculation  of 
the  step  size  if  the  relative  errcr  is  exceeded  for  that 
step.  The  gamma  value  (y  )  is  the  fraction  of  the  old 

re! 

increment  (INC)  step  size  predicted  to  satisfy  the  error 
limitation.  Because  this  method  is  a  fourth  order 
estimation,  the  fraction  required  will  be  the  fourth  root  of 
the  error  equation  (8.18)  and  will  be  defined  as: 
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The  new  step  size  will  be  defined  as: 


INO(new) 


0.8  y 
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INC 
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As  needed,  the  program  will  recalculate  the  step  size  until 
the  size  is  small  enough  to  satisfy  the  error  limitation. 
The  program  for  this  numerical  analysis  is  located  in 
appendix  B. 


The  boundary  conditions  for  this  differential  equation 
are  not  both  initial  value  conditions.  Because  they  contain 
an  initial  and  final  value,  a  'shooting  method"  approach  was 
used  for  the  numerical  method.  An  estimation  for  the 
initial  derivative  of  N  at  R=1  is  predicted.  This 
prediction  is  continually  changed  until  the  final  boundary 
condition  is  obtained. 


8.4  Computer  Simulation. 

Simulations  were  done  using  the  above  method  (see 
appendix  B).  One  simulation  looked  at  the  effects  of 
recombination  while  the  other  looked  at  the  effects  of 
ionization.  The  step  sizes  for  the  simulations  were  reduced 
until  the  step  size  required  for  the  no  generation 
simulation  resulted  in  a  solution  matching  the  analytical 
solution  of  equation  2.49. 

In  the  first  simulation,  the  ionization  frequency  was 

set  to  zero.  Four  runs  of  the  program  were  done  varying  the 

dimensionless  recombination  coefficient  from  0  to  1000  (or 
-  y  3 

irom  0  to  1x10  recombinations  cm  /  sec  in  dimensional 
form).  For  low  recombination  coefficients,  the  plasma 
electrons  extend  to  the  wall.  As  seen  in  Figure  8.1,  higher 
values  of  recombination  cause  the  plasma  to  contract.  As 


figure  8. 1  Radial  Electron  Density  Gradient 
Recombination  Effects 


0.4H 
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observed  in  chapter  six,  the  plasma  width  becomes  smaller  as 
pressure  increases.  This  simulation  illustrates  why  that 
happens.  Additionally,  the  difference  between  strong  and 
weak  regions  can  be  implied  by  the  shape  of  those  curves  (in 
Figure  8.1).  The  curve  clearly  shows  that  the  center  of  the 
plasma  has  the  highest  density  of  electrons.  The  cuter 
region  becomes  less  dense.  Depending  upon  the  cutoff 
density  for  each  observed  region,  an  electron  density- 
profile  may  be  generated  for  plasma  conditions  observed  in 
chapter  six . 

The  second  simulation  held  the  dimensionless 
recombination  at  1000  and  varied  the  ionization  frequency 
from  0  to  30  ionizations  /  sec.  Because  of  the  large  amount 
of  energy  required  to  ionize  a  neutral  atom,  multiple 
collisions  are  required  by  the  electrons  for  ionization  to 
occur.  Therefore,  the  ionization  is  expected  to  be  small 
and  almost  negligible.  As  seen  in  Figure  8.2,  the  electron 
density  near  the  wall  becomes  slightly  larger  when 
ionization  occurs. 

These  simulations  illustrate  that  an  accurate  model 
describing  the  electron  distribution  must  include  the 
generation  term.  These  simulations  do  not  provide  a  very 
accurate  model  as  they  used  constant  values  for 
recombination  and  ionization.  These  values  are  a  function 
of  temperature  arid  pressure,  which  in  turn  depends  upon 
position.  Nevertheless,  the  simulations  do  provide  some 
insight  into  wnat  the  electron  gradient  is  and  why  it  is  for 
various  situations. 


ELECTRON  DENSITY  VERSUS  RADIAL  POSITION 
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Figure  6.? 


Radial  Electron  Density  Gradient  - 
Ionization  Effects 


Radial  Position 


CHAPTER  IX 


CONCLUSIONS 


Several  different  conclusions  can  be  drawn  upon  the 
work  of  this  thesis.  Potential  use  of  this  type  of 
thruster  is  important.  Also,  the  need  for  conducting  a 
microscopic  theoretical  analysis  for  plasma  phenomena  is 
addressed.  The  three  diagnostic  techniques  discussed 
provide  needed  information  to  'h-'-acterize  the  macroscopic 
energy  transfer  within  a  mic  •  nduced  piasma. 

Finally,  computer  numerical  r..o  provide  useful 

information  concerning  plasma  ar.  port  phenomena. 

Because  the  microwave-induced  plasma  e 1 ec t ro t he i mal 
rocket  system  experimentally  displays  similar 
characteristics  to  other  electrothermal  rocket  system?,  it 
has  shown  its  potential  application  for  spacecraft 
propulsion.  Electrothermal  thrusters  lack  high  thrust,  but 
have  high  specific  impulse  which  is  very  useful  for  platform 
station  keeping. 

Investigating  the  microscopic  theory  of  plasma 
particles  is  very  important.  From  the  conservation  of 
particles  and  momentum,  the  shape  of  the  plasma  can  be 
predicted.  From  the  conservation  of  energy,  the  gradients 
of  important  parameters  can  be  predicted,  such  as 
temperature  and  pressure  gradients.  Therefore,  modeling 
transport  properties  must  include  a  microscopic  theoretical 
analysis . 
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The  energy  distribution  experiment  provided  us  with 
the  importance  behind  the  condition  of  the  containment 
wail.  A  dirty  wail  can  absorb  a  large  amount  of  power. 
Reversing  the  gas  flow  showed  a  negligible  difference  in 
the  energy  distribution.  The  power  distribution  appears  to 
stabilize  for  flow  rate  beyond  500  SCCM.  Finally,  initial 
power  absorption  to  the  plasma  is  67%  in  the  TM  mode  and 

80%  in  the  TM  mode. 

0  1  1 

The  plasma  dimensions  experiment  provided  us  with  the 
effects  that  pressure,  flow  rate,  and  power  have  upon  the 
volume  and  shape  of  the  plasma.  The  volume  is  decreased  as 
pressure  increases,  flow  rate  decreases,  or  power  decreases. 
Additionally,  high  flow  rates  tend  to  elastically  elongate 
the  shape  of  the  plasma.  Unlike  nitrogen  plasmas  exhibiting 
an  ellipsoidal  shape,  helium  plasmas  exhibit  a  "dumbbell" 
s hape  . 

The  spectroscopy  experiment  provided  us  with  an 
estimate  for  the  electron  temperature  and  the  effects 
pressure  had  upon  that  temperature.  The  electron 
temperature  is  approximately  13,000  K  in  the  atmospheric 
pressure  region.  Finally,  the  electron  temperature 
decreases  as  pressure  increases. 

The  computer  model  provided  us  with  the  importance  of 
the  recombination  and  ionization  effects.  For  low- 
pressures  when  recombination  and  ionization  coefficients 
are  negligible,  the  wall  recombination  plays  a  very 
important  role  in  the  shape  of  the  plasma.  The  plasma 
reaches  the  wall.  For  high  pressures  when  those 
coefficients  are  important,  wall  recombination  becomes  less 
important.  The  plasma  forms  an  ellipsoidal  shape  and  does 
not  reach  the  wal 1 . 


CHAPTER  X 


RECOMMENDATIONS 


Improvements  can  be  made  on  future  research  in  the 
areas  of  equipment  modification,  experiments,  and 
theoretical  modeling. 

10.1  Equipment  Modifications. 

An  accurate  model  for  the  plasma  region  must  account 
for  the  flow  pattern  (i.e.  the  velocity  profile).  A  major 
area  which  needs  to  be  addressed  is  flow  in  and  out  of  the 
plasma.  For  example,  does  the  flow  by-pass  or  pass  through 
the  plasma  discharge?  If  both  occur,  how  much  of  the  flow 
by-passes  the  discharge?  To  answer  Lhese  questions, 
modification  to  the  collars  c  "  the  plasma  tubes  should  be 
done.  The  collars  should  be  modified  such  that  the  flow- 
velocity  contains  an  angular  component  (swirling  flow). 
Furthermore,  this  modification  should  include  an  adjuster 
such  that  the  axial  helix  interval  can  be  varied. 
Additionally,  a  modification  to  the  collars  should  include 
a  way  for  an  additional  gas  to  enter.  This  additional  gas 
can  be  used  to  conduct  streamline  and  pulse  change 
expe  r  inien  ts . 

Another  important  problem  with  modeling  plasmas  in  a 
micr'ave  resonance  cavity  deals  with  power  absorption  to 
the  walls.  For  an  efficient  thruster  system,  power 
absorbed  by  the  walls  must  be  minimized.  As  seen  in 
chapter  5,  the  condition  of  the  cavity  wall  influences 
power  absorption  to  the  wall.  Two  main  types  of  energy  are 
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absorbed  by  the  wall  -  radiation  and  convection.  To 
minimize  radiation  energy  to  the  wall,  the  inside  wail  of 
the  cavity  should  act  as  a  mirror.  A  mirror  coating,  such 
as  silver,  should  reflect  most  of  the  radiation  energy  to 
the  plasma. 

When  the  investigations  of  analyzing  the  plasma  region 
and  the  thruster  region  result  in  an  accurate  model,  one 
should  consider  combining  the  two  experiments.  For 
example,  a  nozzle  should  be  connected  to  the  cavity  system 
to  investigate  thrust  and  efficiency  measurements.  This 
should  result  in  optimizing  the  location  of  the  plasma 
relative  to  the  nozzle. 

The  quartz  tube  provides  a  large  unknown  temperature 
gradient.  Presently,  we  can  estimate  the  outer  wall 
temperature  of  the  inner  tube  as  being  the  temperature  of 
the  exiting  air  coolant.  An  infrared  temperature  probe  can 
be  used  to  measure  an  intermediate  level  temperature  of  the 
plasma.  This  temperature  can  be  approximated  to  be  that  of 
the  inner  wall  temperature. 

The  wall  effects  on  a  plasma  are  important.  One  such 
effect  is  wall  recombination.  This  wall  recombination 
affects  the  size  and  shape  of  tne  plasma.  To  investigate 
this  effect,  one  should  conduct,  experiments  using  various 
sizes  and  shapes  for  the  confinement  tube. 

Finally,  with  part  of  the  experimental  system  being 
automated,  one  should  consider  automating  the  entire 
system.  The  next  step  in  automating  the  system  includes 
connecting  both  the  temperature  probes  from  the  calorimetry 
experiments  and  the  power  meters  (incident  and  reflective) 
from  the  microwave  power  source  to  the  data  acquistion 
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system.  Once  this  is  completed,  a  computer  program  should 
be  written  to  numerically  analyze  the  data  as  it  is  being 
collected . 

10.2  Proposed  Experiments. 

Several  new  areas  of  experiments  should  be  considered. 
Some  of  them  involve  the  equipment  modifications  as 
mentioned  above.  Dependent  upon  those  modifications,  the 
following  experiments  can  be  conducted.  Flow  pattern 
experiments  can  be  done  to  determine  the  flow  to  and  from 
the  plasma  as  a  function  of  the  environmental  variables 
(i.e.  pressure  and  flow  conditions).  Power  absorption  to 
the  walls  of  the  mirror  coated  cavity  should  be  compared  to 
the  mirrorless  cavity.  Simultaneous  thruster  and  plasma 
experiments  should  be  conducted.  Estimates  of  the  inner 
wall  temperatures  should  be  compared  to  changes  in  the 
environment.  Finally,  simultaneous  experiments  should  be 
done  using  a  completely  automated  system  to  provide  more 
accurate  data  for  constant  environmental  conditions. 

Experiments  which  can  presently  be  done  involve  power 
increases  and  gas  mixtures.  Increasing  the  power  to  above 
one  kilowatt  should  be  done.  This  would  provide  a  more 
realistic  model  characterizing  the  dependence  upon  power. 
Finally,  gas  mixtures  should  be  done  starting  with 
non-reactive  binary  gases.  This  work  should  continue  into 
investigating  tertiary  reactive  gases.  This  evolution  of 
experiments  may  someday  branch  off  into  a  new  type  of  rocket 
thruster  mention  by  Pollard15.  This  new  type  of  thruster  is 
a  hybrid  chemical  /  electric  thruster.  The  reactive  gases 
could  simulate  a  chemical  thruster  with  high  thrust,  while 
the  non-reactive  gases  could  simulate  an  electric  thruster 
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with  high  impulse.  This  is  a  new  concept  which  should 
seriously  be  considered. 

10.3  Theoretical  Modeling. 

Given  the  large  amount  of  data  collected,  work  should 
be  done  in  correlating  the  analyzed  variables.  For 
example,  we  would  need  a  model  to  estimate  the  power 
distribution  and  plasma  dimensions  for  a  given 
environmental  condition  (i.e.  pressure  and  flow  pattern). 

An  initial  step  towards  modeling  the  plasma  region  was 
taken  in  chapter  eight.  The  weakness  of  that  model  stems 
form  its  unrealistic  assumptions  (constant  temperature  and 
pressure  throughout  the  plasma  discharge).  A  more  accurate 
model  should  be  developed  which  accounts  for  the  non-zero 
temperature  gradient.  Additionally,  because  the  electron  is 
not  the  only  important  species  in  the  plasma,  the  gradient 
profile  of  the  other  species  (i.e.  ions  and  neutrals)  should 
be  determined. 

Another  important  model  involves  energy  and  species 
transfer  to  and  from  the  plasma  region.  A  model  should  be 
developed  to  estimate  the  energy  and  mass  transport 
properties  between  the  plasma  and  its  adjoining  fluid. 

10.4  Recommendation  List. 

The  following  summarizes  and  lists  the  above 
recommend at  ions : 

1.  To  modify  the  plasma  tube  collars  to  allow 
experiments  to  be  conducted  for  investigating  flow  patterns 
near  the  plasma  discharge. 
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2.  To  mirror  coat  the  resonance  cavity  for  minimizing 
power  absorption  to  the  wail. 

3.  To  combine  the  cavity  and  thruster  experimental 
systems  for  analyzing  the  optimal  location  of  the  plasma 
with  respect  to  the  nozzle. 

4.  To  use  an  infrared  temperature  probe  for 
estimating  the  inner  wall  temperature. 

5.  To  use  different  size  containment  tubes  for 
investigating  the  wall  effects  upon  the  plasma. 

6.  To  completely  automate  the  experimental  system  so 
that  one  can  conduct  simultaneous  experiments. 

7.  To  conduct  increased  power  experiments  for 
providing  data  with  power  above  one  kilowatt. 

8.  To  conduct  mixture  experiments  for  providing  an 
expanded  data  base  beyond  pure  gas  experiments. 

9.  To  correlate  the  given  and  future  data  for 
predicting  power  distribution  and  plasma  dimensions. 

10.  To  develop  a  theoretical  model  for  determining  the 
species  and  temperature  gradients  within  the  plasma 
discharge . 

11.  To  develop  a  theoretical  model  for  predicting  the 
transport  properties  between  the  plasma  and  the  adjacient 
fluid. 
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APPENDIX  A 
MATHEMATICS 


A.l  Introduction. 

The  following  appendix  highlights  the  vector 
operations,  phasor  transformations  (from  time  domain),  and 
series  solutions  /  orthoganalility  used  in  this  thesis.  Only 
cartesian  and  cylindrical  coordinates  will  be  discussed.  It 
is  noted  that  spherical  coordinates  can  be  used,  but  are  not 
used  in  this  paper. 

3  6 

A. 2  Vector  Operations 

The  following  expressions  and  equations  are  used 
throughout  the  presentation  of  this  thesis: 

Cartesian  Coordinate  Vector: 

E  =  Ex  +  Ey  +  Ez  A.l 

x  y  z 

Cylindrical  Coordinate  Vector: 

E  =  E  r  +  E +  E  z  A .  2 

r  a  z 

Vector  Addition: 

E  +  H  =  ( E  +  H  ) x  +  ( E  +  H  ) y  +  ( E  +  H  )z  A. 3 

xx  y  y  z  z 
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Scalar  Vector  Multiplication: 

—  AAA 

o E  =  HE  x  ♦  <7E  v  +  oE  z 

x  y  z 

Vector  Dot  Product: 

E  •  H  =  EH  +  EH  -cEH 
xx  y  y  z  z 

Vector  Cross  Product: 


E*H  =  (EH  -  EH  )  x  +  (EH-EH  )  y  +  (E  H  -E  H  )z 

y  2  z  y  zxxz  xyyx 


Gradient  of  a  function: 


V  f  j  x  ,  y  ,  z  ) 
vf  ( r  ,  0 , z ) 


df  a  df  a  df  a 

— X  ♦  -3T-y  + 


df  a 
dr  r  + 


dy  J  '  dz 

df  ft  df  a 

r  dd  9  dz  2 


Divergence  of  a  Vector: 

dE 


dE 


dE 


V-E(x, y,z  ) 

V*  E(  r ,0, z  ) 


dx 


d  (  rE  ) 


dy  dz 

dE, 


dE 


dr  +  r  dd  +  dz 


Curi  of  a  Vector: 

f dE 


V*E(x,y,z)  = 


7*E( r,d, z  )  = 


dE 

_ i  _  _ y 

dy  dz 


A 

X  ♦ 


fdE  dE 

X  2 

dz  dx 


A 

y  + 


dE 
_ ) 

dx 


dE 

dy~ 


fdE 

z 

A 

r>  X 

fdE 

r 

dE  1 

Z 

§  «.  i 

d ( rEg ) 

dE  ') 
r 

rdd 

dz 

j 

r  t 

dz 

h. 

dr 

> 

r 

dr 

de 

J 

A.4 

A .  5 

A .  6 

A .  7 

A. 8 

A.  9 

A.  10 

A.  11 

A  .  12 
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Laplacian  of  a  function: 


V2f (x, y , z ) 


7  f  (  r  ,  0  ,  z  ) 


=  ?*<*f) 


a2f 

+  §h  + 

a2f 

dx'a 

dy2 

dz2 

d'f  ♦  - 

d2t 

+  — 

r 

2  2 

66 

dzZ 

A.  1  3 


A.  14 


A. 3  Phasor  Transformation, 


Phasors  are  defined  as  transformation  of  a  differential 
or  integral  of  a  time  dependant  function  into  an  algebraic 
expression.  To  accomplish  this,  one  must  rewrite  the  time 
domain  function  as  the  real  component  of  a  complex  phase 
domain  function. 


Let : 


f ( t )  =  cosiut ) 

Rewritten  as  real  component  of  complex  phase  domain 
f unct ion : 

f(t)  =  Re^cos(ot)  ♦  j  sin(wt)J 

Differential 


Di  f  f erent ial 

d  f  (  t  ) 
dt 

d  f  (  t  ) 
dt 

d  f(t) 


equat ion : 


R  fd  cos  hit  )  d  [  ,j  s i n ( ot  )  ] 
eL  dt  +  dt 


Re 


-u 


sinlwt )  +  jo  cos(ot) 


A  .  1  5 


A.  16 


A.  16 

A.  17 


dt 


-u  s in ! ut  ) 


A.  18 


Define  Differential  Phasor: 


luo 


Substitution : 

A.  28 

A.  28 

A.  20 

A.  4  Series  Solutions  and  Orthoganal ity . 


[f ( t  )  dt 

J  n  t )  d  t 

!  f  <  t  >  d  t 


Re  cos(ut)  +  j  sin(wt) 


Re  T  cos(ut)  ♦ 

l  Ju 


sinfwt ) 


s  in  ( <Jt 
w 


A  series  solution  method  is  used  to  calculate  the 
solution  of  special  kinds  of  higher  order  differential 
equations.  For  this  thesis,  only  the  Bessel  differential 
equation  is  used.  Or thoganal l ty  is  used  for  calculation  the 
coefficients  ot  the  solution.  Bessel  functions  and  certain 
types  of  trigonometric  series  are  orthoganal. 


Bessel's  Differential  Equation: 


2  d  y  dy  2  2 

x  — §-  ♦  x  »  ( X  -  V,  1  y 

ax 


A.  31 


Series  Solution  (Frobenius  Method)  to  equation  A. 31 


y  i  x  1  =  ;  a  x 

m 

m«  0 


I 


ra+V 


A.  32 


Obvious  Solution  for  V  : 


V  = 


A.  33 


Define  Bessel  Function  of  First  Kind 


vx> 


■’I 

ra*  o 


t-i)"1  *2ni 
22*''1*’  m  !  (  iftm  )  ! 


Substitution  for  V  : 


y  (  x  )  -  V 


/ 

n  ■  0 


C  J  (ax)  + 
n  U)  n 
1 


D  J  ( a  x  > 
n  -v,  n  J 


Solutions  for  C  and  D  can  be  calculated  from  the 

n  n 

properties  of  or t hoganal i t y . 

Define  Inner  Product  of  Two  Functions: 


<  f  (  x  )  ,  g  (  x  )  > 


g(  X  ) 


dx 


Inner  Product  of  Or t hog  anal  Functions: 

<  f ( x ) , g ( x ) >  =  0  for  f ( x  )  *  g ( x  ) 

Rewriting  y( x  )  : 

y  (  x  )  =  C  f  (  x  )  +  C  g  (  x  )  +  °  ■>  ’ 

To  solve  for  C  ,  one  must  multiply  each  side  by  fix)  and  a 


A.  34 


A  .  3  5 


A  .  36 


A  •  3  i 


A  .  38 
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weight  function,  w ( x  i  ,  and  integrate  over  the  domain  of  the 
f unc  t ion . 


L 


w(  X  )  f  (  X  )  y  (  x  )dx  =  Jc^w  (  x  )  f  (  x  i  f  (  x  )dx  +  Jc^w  i  x  )  f  (  x  ) g  i  x  ) dx  ♦  A.  39 


By  orthoganal ity  of  the  series: 


j'w(x)  f  (  x )  y(x)  =  Jc^wtx)  f(x)  fix)  dx  +  0 


A.  4  0 


Rearranging  the  terms 


I  W 


(  x  )  fix)  ylx  i  dx 


I  “ix 


)  f ( x )  f ( x )  dx 


A.  4  1 


iome  useful  integrals  of  Bessei  functions 


37 


fx"  J  (x)  dx  =  x n  J  (  x)  +  C 

J  n-1  r\  1 


A. 42 


jJn*l<X)  dX  = 


J  ,  (x)  dx  -  2J  (x) 

n-1  n 


A. 43 


\nJ  (x)  dx  =  -x  "  J  (x)  +  C 

ii  ♦!  n  1 


A.  44 


fx  J2(x)  dx  =  ~~  J2  ,  (  R ) 


2  n  1 1 


A.  45 


o 

39 

A, 5  Useful  Vector  Properties 

The  following  are  several  useful  vector  relationships 
which  become  very  helpful  when  solving  vector  related 
equations  such  as  Maxwell's  equations  for  electromagnetics. 

A  (B*C)  =  (AxB)  C  =  B  {C*A)  A. 46 
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A  MB  x  C)  =  |A  1  C)  B  -  (A  '  B)  C  A. 47 

(A  *  B)  •  iC  *  D)  =  !  A  •  C)  (B  ■  D)  -  (A  •  D)  (B  •  C>  A. 48 

7  i f  g )  =  f  Vg  +  g  Vf  A. 49 

?  •  <  f  A)  =  f  7-a  +  A  ■  Vf  A.  50 

7  x  if  A)  =  f  7xa  +  Vf  x  a  A. 51 

72~  =  V  (7- A)  -  V  x  (VxA)  A.  52 

V  x  Vf  =  V  •  (Vxa)  =  0  A. 53 


APPENDIX  B 


*  PROGRAM 

*  NAME 

*  DATE 
•DESCRIPTION 


M:  Runge  Kutta  Solution  for  Electron  Diffusion  * 

E:  Scott  S.  Haraburda  * 

E:  1  January  1990  * 

N:  This  program  is  designed  to  estimate  electron  * 
concentration  gradient  within  the  plasma  region  * 
using  the  fifth  order  Runge  Kutta  numerical  method* 
using  variable  step  size  control.  This  control  is* 
accomplished  using  an  absolute  error  test.  * 

* 

VARIABLES : HSTART  -  The  Initial  Step  Size  for  the  Program.* 

HMIN  -  The  Minimum  Step  Size  for  the  Program.* 

HMAX  -  The  Maximum  Step  Size  for  the  Program.* 


HMAX 

RTOL 

RERR 

G 

R 

RQ 

RR 

NO 

NR 

DT 

ROLD 

KNEW 

NOLB 

NNEW 

GAMMA 

I  FLAG 

RK1-6 

NK1-6 

EST 

A 1  -5 

Rl-6 

I 

IR 


-  The  Relative  Tolerance  Level. 

-  The  Relative  Error  for  a  Given  Step. 

-  The  Initial  N  Value  in  the  Interval. 


-  The  Initial  N  Value  in  th< 

-  The  Latter  N  Value  in  the 

-  The  Value  of  R  at  Point  Q. 

-  The  Value  of  R  at  Point  R. 

-  The  Value  of  N  at  Point  0. 

-  The  Value  of  N  at  Point  R. 


Interval . 


at  Point 
at  Point 
at  Point 
Size  for 


-  The  Value  of  N  at  Point  R. 

-  The  Given  Step  Size  for  the  Program 

-  The  Old  Value  of  R. 

-  The  New  Value  of  R. 

-  The  Old  Value  of  N. 

-  The  New  Value  of  N. 

-  The  Calculated  Fractional  Step  Size 

-  The  Termination  and  Write  Flag. 

-  The  Runge  Kutta  Functions. 

-  The  Runge  Kutta  Functions. 

-  The  Estimated  Error. 

-  The  Runge  Kutta  Coefficients  for  New 

-  The  Runge  Kutta  Coefficients  for  EST. 

-  The  Loop  Integer  for  Main  Program. 

-  The  Collision  Rate  Frequency  Value 


Step  Size 
Flag. 


RECOMB  -  The  Recombination  Rate  Frequency  Value* 
REST  -  The  Estimated  Value  for  R.  * 
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★ 

NEST 

-  The 

Estimated  Value 

for  N . 

* 

★ 

NT 

-  The 

intermediate  New 

Value 

of  N. 

* 

* 

RT 

-  The 

Intermediate  New 

Value 

of  R. 

* 

*  * 
************************************************************** 

I  inTEGER  I,  I  FLAG 

REAL  HSTART , HMIN , HMAX , RTOL , RERR ,  Q ,  R , DT , RECOMB , I R , REST , 

+  NERR , NEST 

DOUBLE  PRECISION  RQ ,  RR ,  NQ ,  NR 

OPEN  (9,  FILE  =  ’ ARMY 1 ’ ,  STATUS  =  ’NEW’) 

1  PRINT  *,  'Input  the  Recombination  Coefficient  [RECOMB):’ 

READ  *,  RECOMB 

PRINT  *,  'Input  the  Collision  Frequency  [ I R ) : * 

READ  *,  IR 

PRINT  *,  ’Input  the  Initial  Derivative  Value  (NQ):’ 

READ  * ,  NQ 

PRINT  * ,  ’Do  You  Want  the  Results  Saved  [1-yes,  u-no|?’ 

READ  *,  I  FLAG 

HSTART  =  0.01 

HMIN  =  IE-6 

RTOL  =  IE-6 

HMAX  =  0.01 

R  =  1  . 

RQ  =  0. 

DT  =  0.01 
PRINT  75 

PRINT  100,  R,  NQ,  RQ ,  R,  R,  HSTART 
<5=1. 

IF  (I  FLAG  .EQ.  1)  WRITE  (UNIT  =  9,  FMT  =  *)  Q,  RQ 
DO  5  I  =  1,99 

R  =  1  .  -  (  I  -  1  )  *  DT 

Q  =  1  .  -  I  *  DT 

CALL  DESOLV  ( Q , R , NQ , NR , RQ , RR , RTOL , HSTART , HM IN , HMAX , 

+  IFLAG, NERR, RERR, RECOMB, IR) 

IF  ( IFLAG  . EQ.  3 )  GOTO  10 

PRINT  100,  Q,  NR,  RR,  NERR,  RERR,  HSTART 

IF  (IFLAG  .EQ.  1)  WRITE  (UNIT  =  9,  FMT  =  * )  Q ,  RR 

NQ  =  NR 

RQ  -  RR 

5  CONTINUE 

IF  ( IFLAG  .EQ.  0 )  GOTO  1 
GOTO  50 

10  PRINT  *, ’Algorithm  failed.  Last  point  was  I R , NOLD , ROLD ) 

PRINT  *,  Q,  NR,  RR 
50  CONTINUE 

75  FORMAT  ( / / , T5 , ’ r ’  , T 1 4 ,  ’ N ( r ) ’ , T28 , ’ R( r  )  ’ , T3 9 ,  ’ N  Rel  Err’, 

+  T 5 3 ,  ’ R  Rel  Err ’  , T66  ,  ’ HSTART ’  , / , 7 3 ( 1H-  )  ) 

100  FORMAT  ( T2 , F5 . 2 ,T9 , E12 . 5 ,T23 ,E12 . 5 ,T37 , E12 . 5 , T51 ,E12 . 5 , 

+  T65.F7.4) 

CLOSE  (9) 

END 
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SUBROUTINE  RKF  (NOLD,  NNEW,  ROLD,  RNEW ,  NEST,  REST,  INC, 
+  RECOMB,  IR,  R) 


*  * 

*  The  Runge  Kutta  fourth  /  fifth  order  method  to  estimate  * 

*  rhe  next  vaiue  of  N i r )  4  k(r)  and  their  estimated  error.  * 

*  * 


INTEGER  I,  I  FLAG 

REAL  NEST, REST, INC ,NK1 , NK2 ,NK3 , NK4 , NK5 , NK6 , Q , R , DT , IR, 
+  RK1  ,  RK2  ,  RK3  ,  RK4  ,  RK.5  ,  RK6  ,  RECOMB ,  RTOL  ,  RERR ,  NERR ,  RN 

DOUBLE  PRECISION  ROLD , RNEW , NOLD , NNEW , NT , RT 

*********************************************************** 


X 


* 


*  Define  the  Runge  Kutta  coefficients.  * 

*  * 


B2  1 

= 

i  .  / 

4  . 

B22 

= 

3  .  / 

8  . 

B23 

- 

1  2  . 

/  13. 

B31 

3.  / 

32  . 

B3  2 

9.  / 

32  . 

B41 

= 

1932 

.  /  2197 

B42 

r 

-7200.  /  219 

B43 

= 

7296 

.  /  219  7 

B51 

= 

439. 

/  216. 

B52 

= 

-8. 

B53 

= 

3680 

.  /  513  . 

B54 

- 

-845 

.  /  4104 

B61 

-8. 

i  2  7. 

B62 

2  . 

B63 

= 

-3544  /  2565 

B64 

= 

1859 

.  /  4104 

B65 

- 

-11 . 

/  40. 

Al 

= 

25  .  / 

216. 

A  3 

1408. 

/  2565. 

A  4 

= 

2197. 

/  4104. 

A  5 

= 

-1  .  / 

5  , 

R 1 

= 

1  .  / 

360. 

R3 

= 

-128. 

/  4275  . 

R4 

-2197 

.  /  7524 

R5 

1  .  / 

50  . 

R6 

= 

2.  / 

5  5  , 

NT 

= 

NOLD 

RT 

= 

ROLD 

112 

************************************************************* 
*  * 

*  Begin  the  Runge  Kutta  numerical  analysis.  * 

*  * 
************************************************%********£*** 

RN  =  R 

N'Kl  =  -NT/RN  +  RECOMB*RT*RT  -  Ift*RT 
RK1  =  NT 

NT  =  NOLD  +  INC  *  B21  *  NK1 
RT  =  ROLB  +  INC  *  B21  *  RK1 
RN  =  R  +  INC  *  B21 

NK2  =  -NT/RN  +  RECOMB* RT*RT  -  IR*RT 
RR2  =  NT 

NT  =  NOLB  +  INC  *  (B31*NK1  +  B32*NK2 > 

RT  =  ROLB  +  INC  *  (B31*RK1  +  B32*RK2  ) 

RN  =  R  +  INC  *  B22 

NK3  =  -NT/RN  +  RECOMB*RT*RT  -  I R*RT 
RK3  =  NT 

NT  =  NOLB  +  INC  *  (B41*NK1  +  B42*NK2  +  B43*NK3) 

RT  =  ROLB  ♦  INC  *  (B41*RK1  +  B42*RK2  +  B43*RK3) 

RN  =  R  +  INC  *  B23 

NK4  =  -NT/RN  +  RECOMB*RT*RT  -  IR*RT 
RR4  =  NT 

NT  =  NOLB  +  INC  *  ( B5 1 *NK1 +B52*NK2+B5 3*NK3 +B54*NK4 ) 

RT  =  ROLB  +  INC  *  ( B5 1 *RK1 +952*RK2+B5 3*RK3 +B54 *RK4 ) 

RN  =  R  +  INC 

NK5  =  -NT/RN  +  RRCOMB*RT*RT  -  I R*RT 
RK5  =  NT 

NT  =  NOLB  +  INC*( B61*NK1+B62*NK2+B63*NK3*B64*NK4+B65*NK5  ) 
RT  =  ROLB+INC*(  B61*RK1  +  B62 *RK2  +  B63*RK3+B64*RK4+B65*RK5  ) 
NNEW  =  NOLD  +  INC  *  <A1*NK1  +  A3*NK3  ♦  A4*NK4  +  Ao*NK5 ) 
RNEW  =  RuLD  +  INC  *  <A1*RK1  +  A3*RK3  +  A4*RK4  +  Ao*RK5 ) 
RN  =  R  +  INC  /  2. 

NK6  =  -NT/RN  +  RECOMB*RT*RT  -  IR*RT 
Rh.6  =  NT 

NEST  =  R 1  *  N  K 1  +  R3*NK3  +  R4*NK4  +  R5*NK5  +  Rd*NK6 
REST  =  R1*RK1  +  R3*RK3  +  R4*RK4  ♦  R5*RK5  +  Rb*RK6 
RETURN 
END 
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SUBROUTINE  GAMCAL  ( Q , R , NOLD , ROLD , NEST , REST , RTOL , GAMMA ) 


*  The  method  to  evaluate  the  necessary  step  size  for  the  * 

*  given  tolerance  levels  and  precalculated  estimated  error.  * 

*  The  mixed  control  is  evaluated  here  using  both  the  relat.  * 

*  and  absolute  tolerances.  * 


DOUBLE  PRECISION  NOLD, ROLD 
REAL  T 1 , T2 , ABSEST , NEST , REST , 0 , R , RTOL 
IF  (NEST  .NE.  0.)  THEN 
T 1  =  ABS  (NOLD  /  NEST) 

ELSE 

T 1  =  1 E30 
END  IF 

IF  (REST  .NE.  0.)  THEN 
T 2  =  ABS  (ROLD  /  REST) 

ELSE 


T2  =  1 E30 
END  I  F 

ABSEST  =  MIN  (T1,T2> 


*  Calculate  the  gamma  value  for  the  variable  step  algorithm.* 

*  * 


************************************************************* 

IF  (ABSEST  .EQ.  0. )  GOTO  1 

GAMMA  =  (RTOL  *  ABSEST  /  ( R  -  Q >  )  **  0.25 

RETURN 

1  Ti  =  MAX  ( NEST, REST, RTOL/ 10  ) 

GAMMA  =  (RTOL  *  RTOL  /  (Tl  *  (R  -  0)))  **  0.25 

RETURN- 

END 
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SUBROUTINE  DESOLV  ( Q , R , NQ , NR , RQ , RR , RTOL , HSTART , HMIN , 

+  HMAX  ,  I  FLAG, NERR.RERR, RECOMB,  IR) 

************************************************************ 
*  * 

*  The  main  subroutine  combining  the  Runge  Kutta  algorithm  * 

*  along  with  the  variable  step  size  calculation  for  a  giv 

*  interval  [Q,R].  Minimum  and  Maximum  step  sizes  are  nee 

*  This  subroutine  is  terminated  if  the  calculated  step  size* 

*  is  less  than  the  minimum  step  size  provided.  * 

*  * 


INTEGER  I  FLAG 

REAL  Q , R , HSTART , HMIN , HMAX , RTOL , RERR , NERR , NEXT , REST , 

♦  HOLD, BT,T, NEST, RECOMB, IR 

DOUBLE  PRECISION  NQ , NR , RQ , RR , ROLL , RNEW , NNEW , NOLB 
HOLD  =  HSTART 
T  =  R 
NOLD  =  NO 
ROLD  =  RQ 

************************************************************ 

*  * 

*  Using  given  parameters,  the  new  value  and  required  step  * 

*  size  is  calculated.  * 

*  * 
************************************************************ 

1  CALL  RKF  ( NOLD, NNEW, ROLD. RNEW, NEST, REST, HOLD, RECOMB, IR,T) 
CALL  GAMCAL  ( Q , T , NOLD , ROLD , NEST , REST , RTOL , GAMMA ) 

HNEW  =  (0.8)  *  GAMMA  *  HOLD 
IF  (GAMMA  .GE.  1 . )  GOTO  2 

IF  (HNEW  .LT.  HOLD/ 10.)  HNEW  =  HOLD  /  10. 

IF  (HNEW  .LT.HMIN)  GOTO  3 
HOLD  =  HNEW 
GOTO  1 

*********************************************************** 

*  * 

*  Calculated  new  value  for  R  is  acceptable,  and  next  * 

*  parameters  are  established.  * 

*  * 

*********************************************************** 

2  IF  (HNEW  .LT.  HMIN)  GOTO  3 

IF  (HNEW  . GT.  5 . *H0LB )  HNEW  =  5.  *  HOLD 

IF  (HNEW  .GT.  HMAX)  HNEW  =  HMAX 

IF  (T-HOLD  . LE.  Q)  GOTO  4 

T  =  T  -  HOLD 

HOLD  =  HNEW 

NOLD  =  NNEW 

ROLD  =  RNEW 

GOTO  1 


a  m 


Algorithm  failure  flag  and  last  values  are  set. 


*  * 

3  I  FLAG  =  3 

R  =  T 
NR  =  SOLD 
RR  =  ROLD 
RETURN 


Computation  of  final  step  in  interval,  along  with  the 
relative  and  absolute  errors  (of  that  final  step). 


4  CONTINUE 

HSIART  =  HNEW 
HOLD  =  T-Q 

CALL  RKF  ( NOLD , NNEW , ROLD , RNEW , NEST , REST , HOLD , RECOMB , IR , T ) 
NR  =  NNEW 
RR  =  RNEW 

IF  (RR  .  EQ.  0.)  THEN 
RERR  =  1  . 

ELSE 

RERR  =  ABS ( REST  /  RR ) 

ENDIF 

IF  (NR  .EQ.  0. )  THEN 
NERR  =  1 . 

ELSE 

NERR  =  ABS (NEST  /  NR) 

ENDIF 

RETURN 

END 


